Theoretical Modeling and Parameter Analysis of Micro-Pulsed Plasma Thruster

In this paper, a new numerical modeling method that combines a modified electromechanical model with a Teflon ablation model is proposed to simulate the working process of micro-pulsed plasma thruster (μPPT). The ablation mass accumulation during the discharge process and the non-Fourier effect during the ablation process are considered to improve simulation accuracy. Simulation results are in good agreement with experimental results. Furthermore, the influence of both the electrical and structural parts on the performance of μPPT was investigated by applying the combined model. It was proven that large capacitance of the capacitor and large gap-to-plate-width ratios can effectively improve the comprehensive performance of the thruster concluding the specific impulse, impulse bit, and efficiency. The results not only validate the theoretical model, but also guide the design and operation optimization of μPPT.


Introduction
Pulsed plasma thrusters (PPTs) are considered an attractive propulsion option for stationkeeping and drag makeup purposes for mass-and power-limited satellites that require micro-Newton second to milli-Newton second impulse bits [1][2][3].The low power requirements, simple design, robustness, and high specific impulse (>1000 s) are some of the key factors contributing to the application of PPT [4].With plenty of illustrations of the maturity of this technology, it has also been shown that PPT can be easily scaled down in power and size.A micro-PPT (µPPT), which is a miniature version of the traditional PPT, has been designed for delivery of very small impulse bits.The µPPT can deliver an impulse bit in the 10 µN•s range to provide attitude control and stationkeeping for microsatellites.In this thruster, the discharge across the propellant surface ablates a portion of the propellant, ionizes it, and then accelerates it, generating the thrust in a predominantly electromagnetic way [1,4,5].
There is an increasing need for precision from the control of micro/nano satellites, which can be achieved by optimizing µPPT processes.Since experiments are cost-intensive, numerical methods are often used to simulate the working processes of the µPPT, so as to explore the effects of relevant parameters on improving performance [6][7][8].
The modeling of µPPT has many aspects in common with the simulation of PPT, which has recently been extensively reported.Among these models, the electromechanical model is one of the most frequently used in the simulation of PPT.This model approximates the PPT system as an electromechanical device with an electric circuit interacting with a dynamic system.The electric circuit is theoretically idealized in both models as a Lenz-Capacitor-Resistance (LCR) circuit with discrete, albeit moveable, elements, as shown in Figure 1.With these assumptions, the electromechanical model is relatively mature and has the advantages of small computation, short computing period, and simple modeling steps compared with the slug model [9] and magnetohydrodynamic (MHD) model [10].Moreover, the simulation results were in good agreement with the experimental results, and the influences of parameters on the performance of µPPT can be investigated through this model.The electromechanical model, which has gained wide application in modeling, was first described by Hart in his fundamental 1962 study [11].Jahn [9] developed a one-dimensional model that can predict the macroscopic properties of PPT with simple experimental results.However, it was inconsistent with the actual situation of increased ablation mass and could not reflect the ablation and plasma flow processes, as the ablation mass was assumed to be a constant value and the ablated working fluid was assumed to be directly attached rather than flowing to the current sheet during the discharge process.Subsequently, scholars have made many improvements to the electromechanical model.Vondra et al. [12] developed a modified model that took the effect of aerodynamic force on the acceleration process of the working fluid into account and applied it to the optimized design of the Lincoln Experimental Satellite 6 (LES-6) PPT.Wei [13] proposed a diffusion model that can reflect variations in the thickness of the current sheet with the discharge current.In these papers, the models were mainly modeled with empirical parameters or experimental results and could predict macroscopic performance rapidly, but they failed to accurately reflect the ablation process.Moreover, they assumed that the ablation mass as a constant value was no longer negligible in the simulation of µPPT due to the low magnitude in volume, power, and weight.It is therefore necessary to establish a theoretical model that contains the ablation process and relates the accumulation of the ablation mass to the mechanism process of impulse formation [14][15][16].
With relatively low density, steady physical and chemical properties, low conductivity, and high specific heat, Teflon has become an ablation-resistant propellant, and is applied in PPT [17][18][19] and other aerospace applications [20].There have been a number of investigations into the ablation process of Teflon, and some models have been presented to help study the ablation problem [21].To survey the transient ablation of Teflon, heat conduction and phase transition in intense radiation and convective environments were investigated by Arai [22].The mechanisms of surface recession and the evolution of the internal temperature and thickness of the gel layer were investigated with the transient one-dimensional two-layer ablation model of Teflon.Arai's work took into account the optical transmittance of the amorphous and crystalline zones of Teflon, and studied the effect of the internal absorption of radiation in a numerical simulation.On the basis of Arai's work, Gatsonis et al. [23] modified the ablation model to simulate PPT operations and developed a volume fraction method to capture the moving interface between the semi-crystalline and gel phases; the simulation results were in good agreement with the experimental results.However, in situations that include extremely high temperature gradients, extremely large heat fluxes, and extremely short transient durations, the propagation speeds are finite and the mode of heat conduction is propagative and nondiffusive.As a result, the effect of non-Fourier heat conduction on the ablation process of Teflon cannot be neglected [24].In particular, the non-Fourier effect would be more obvious under the operation of µPPT due to the lower magnitude.But Gatsonis et al. did not take the non-Fourier effect into account.To analyze the physical process in the Teflon cavity of a coaxial PPT, Keidar [25] also built an ablation model to calculate the electron temperature, electron and neutral densities, and ablation surface temperature with electron temperature and electron and neutral densities, considering the plasma energy balance, Joule heating, thermal conductivity, and mass balance.However, Keidar also paid little attention to the non-Fourier effect and some performance parameters could not be calculated in his paper.Additionally, there is little focus on the non-Fourier effect in similar attempts like those of Cho [26], Huang [17], and Tahara [27].It is therefore significant to take the non-Fourier effect into account in the ablation process.This paper presents a numerical modeling method that combines a modified electromechanical model with a Teflon ablation model.In our method, a modified ablation model is established to calculate the ablation mass accumulation of work fluid, taking into account the reflectivity of material, heat conduction, phase transition, and non-Fourier effect.After this, the cumulative ablation mass is assumed to be an input parameter of the modified electromechanical model to simulate the operation of µPPT.Implementing the combined modeling method, the precision is improved while the calculation increases only very little.Meanwhile, by utilizing the model, not only are the macroscopic performance parameters predicted, but the specific ablation process and the temperature distribution of the working fluid can be given [14,28].
The remainder of the paper is organized as follows.Section 2 presents the theoretical modeling process of the ablation model and modified electromechanical model.Section 3 compares the numerical results with the experimental results in order to verify the effectiveness of the modeling method.Section 4 explores the influences of different parameters on the performance of µPPT.Finally, some conclusive remarks are given in Section 5.

Theoretical Modeling Method
To put more attention on the modeling process, we chose a micro parallel-plate pulsed plasma thruster as the modeling object.This µPPT is a solid-propellant plasma accelerator that comes in a parallel-plate and coaxial configuration, as shown in Figure 1.It consists of two electrodes that are connected to a capacitor and uses solid Teflon as a propellant.A spark plug generates some plasma that provides an initial conductive path between the electrodes.The main capacitor is discharged, generating an arc that ablates and ionizes the Teflon propellant.The generated plasma is accelerated through electromagnetic and gas dynamic forces to produce thrust [15,26,29].

Ablation Model
In the working process of µPPT, the spark plug is used to induce the discharge of the energy storage capacitor.Then, under the discharge arc, the energy gradually accumulates on the ablation surface of the Teflon solid layer and the temperature of the ablation surface increases.Given that the phase-transition temperature of Teflon (T p ) is 600 K, the heat conduction and phase transition processes are described in the following three ablation stages, as shown in Figure 2, according to the different phase-transformation states.
In the first ablation stage, the surface temperature is below T p and phase transition does not exist in the internal and boundary positions.A single-layer thermal model is applied and the heat conduction equation is solved to obtain the temperature distribution of the single Teflon solid layer.When the surface temperature is raised to T p , the phase transition phenomenon occurs, and this is named the second ablation stage.In this stage, a two-layer thermal model is adopted to investigate the transient heat conduction and phase transition in the target, which can be divided into two layers, the solid layer and the gel layer of the Teflon.In particular, a significant melting phase interface forms between the solid and gel layers of the target.The initial length of the Teflon solid layer is labeled L in Figure 1.The location of the melting phase interface is labeled L m , while the locations of the exposed ablation interfaces are labeled L, L 1 , and L 2 .Additionally, the moving velocity of the melting phase interface and exposed ablation interface are v m and v, respectively.The temperature distribution of the two ablation layers can be obtained by solving the heat conduction equations with different boundary conditions, which will be presented below.When the discharge is over, the heat will disappear.Following this, the working substance is in a state of natural cooling.At this third ablation stage, if the temperature is still high, the Teflon solid will continue to be ablated while the ablated plasma cannot get electromagnetic acceleration.This process is also known as delayed ablation [20,21,30].

First Ablation Stage
In the first ablation stage, the Teflon solid layer is heated up, but the temperature is still below the phase transition temperature (T p ).This means that the Teflon solid is in crystalline phase without a phase transition, as shown in Figure 2a.The energy used for the ablation process comes from the arc of the capacitor discharge; assume that the ratio of the single pulse discharge energy absorbed by the surface of the Teflon solid is γ [14,19,20,31].Therefore, the heat flux acting on the surface of the working fluid is expressed as: where V(t) is the capacitor discharge voltage, I(t) is the circuit discharge current, h stands for the plate gap, and w stands for the plate width.The temperature in the Teflon solid layer can be described by non-Fourier heat conduction [20]: where k s is the thermal conductivity of the Teflon under a solid state, ρ s is the density under a solid state, C s is the specific heat capacity of the Teflon solid, V a is the acoustic velocity, and τ 0s is the thermal relaxation time under a solid state.Physically, the thermal relaxation time τ 0 is regarded as the macroscopic parameter of a series of microscopic interactions in the target.It describes the mean time effect of all these microscopic relaxation processes, such as electron-electron, electron-phonon, and phonon-phonon scattering.
An adiabatic boundary condition is applied at the rear surface of the Teflon solid: A heat flux boundary condition is applied at the outer exposed surface of the Teflon solid layer and radiation heat transfer is taken into account: where σ is the Boltzmann constant, ε is the absorption coefficient of the working mass ablation surface, and T s (L,t) is the surface temperature of the Teflon solid [14].

Second Ablation Stage
With the arc continuing to heat the Teflon, the temperature of the ablation surface gradually rises.When the surface temperature surpasses the phase transition temperature T p at time t m , a phase transition process from solid to gel occurs, which we call the second ablation stage.A recessive melting phase interface forms and the Teflon target can be divided into a solid layer and a gel layer by the phase interface, as shown in Figure 2b.With the analysis of Zhang's study, the initial value of the location of the exposed ablation interface equals the initial length of the target, i.e., L = L 1 [21].
As far as the Teflon solid layer is concerned, the temperature is still determined by Equations ( 2) and (3) [21].Meanwhile, the temperature of the Teflon gel layer is described as follows: ) where k m is the thermal conductivity of the Teflon under a gel state, ρ m is the density under a gel state, C m is the specific heat capacity of the Teflon gel, V a is the acoustic velocity, τ 0m is the thermal relaxation time under the gel state, and Q m is the energy released per unit time and volume (W/m 3 ) in the depolymerization process.Q m is given as Equation ( 8): where A p is the pre-exponential factor, E a is the activation energy, R is the universal gas constant, and E p is the specific depolymerization energy [21].
The boundary condition of the ablation surface of the gel layer is given as follows: Furthermore, the Teflon decomposes, with gaseous monomer vapors coming from the exposed surface to the vacuum circumstance.The mass flux from the exposed Teflon surface can be calculated by integrating the depolymerization rate over the entire amorphous domain via [21]: )dx (10) where L m is the location of the phase transition plane at the current time and given as the following equation: where v m is the moving velocity of the phase interface, H 0 is the latent heat of the solid-to-gel phase transition, k m is the thermal conductivity of the Teflon under a melt state, and ρ 0 is the mean density.The recession velocity is expressed as Equation ( 13): where ρ r is the reference density of the Teflon.
In the single pulse time, the moving distance of the ablative surface is defined as Equation ( 14): where t is the time of a single pulse.

Third Ablation Stage
Following this, the working substance in this stage is as shown in Figure 2c.At this stage, there is no heat resource, and if the temperature is still high, the Teflon solid will continue to be ablated while the ablated plasma cannot get electromagnetic acceleration.This process is also known as delayed ablation.
If the temperature of the Teflon is higher than the melting temperature, the internal temperature is described by Equation ( 15) and the boundary condition of the ablative surface is defined by Equation (16).Furthermore, Equations ( 17) and ( 18) correspond to the temperature distribution and the ablative boundary conditions at a lower temperature. ) where L 2 is given as Equation ( 19):

Modified Electromechanical µPPT Model
To simplify the numerical analysis, the following assumptions were adopted: 1.
The parallel plate µPPT thruster, shown in Figure 1, with the moving discharge is described electrically as an LCR circuit.The electric circuit is theoretically idealized in both models as an LCR circuit with discrete, albeit moveable, elements [28].

2.
This model approximates the µPPT system as an electromechanical device with an electric circuit interacting with a dynamic system [26].

3.
The mass of Teflon ablation accumulates with the discharge process in the working process of µPPT, and this is different from the hypothesis of the conventional model, hence higher model accuracy is achieved.Moreover, the ablative mass, m(t), can be obtained by integrating the mass flow in Equation ( 10). 4.
Only the gas dynamic force and the Lorenz force are considered in the force analysis of the current sheet.

5.
The plasma produced by the ablation is in a state of complete ionization.

6.
The ratio of the single pulse discharge energy absorbed by the surface of the Teflon solid, γ, is assumed to be a constant value.7.
The thruster is assumed to be in an ideal and stable working state.
The physical law used to describe the dynamics of the circuit is Kirchoff's voltage law, which is itself derived from Maxwell's charge continuity equation for the circuit.Therefore, the discharge process of the loop can be expressed by the following equations: where R t (t) is the sum of the loop resistance; L T (t) is the sum of loop self-inductance; V 0 is the initial voltage of the capacitor; C is the capacitance value; and L c , L e , and L pe are, respectively, the capacitor self-inductance, the wire and lead self-inductance, and the plate electrode geometry self-inductance.The resistance terms R c , R e , R ps , and R p refer, respectively, to the capacitor, the wires and leads, the plate electrodes, and the plasma [28].
The motion equation of the current sheet can be rewritten as Equation (23): .
where x(t) is the vertical distance between the current sheet and the working mass ablation surface, .
x(t) is the kinematic speed of the current sheet, m(t) is the cumulative mass of the working mass ablation at a given moment, m(t 0 ) is the amount of ablative mass accumulated at the last time step, .m(t) is the ablative mass at every time step, F L is the Lorenz force, and F A is the gas dynamic force.
As reported in [14], F L and F A can be calculated as follows: where h is the distance between two electrodes, w is the width of the electrodes, d is the thickness of the electrodes, n e is the volumic electron number, σ is the Boltzmann constant, T e is the electron temperature, and µ 0 is the permeability of the vacuum.
The plasma caused resistance is given as [14]: e µ 0 ln 1.24 × 10 7 ( T 3 e n e ) Based on the above governing equations, one modified electromechanical µPPT model is developed as follows: The main parameters related to the performance of the thruster are the specific impulse, the impulse bit, the efficiency, and the ratio of thrust to power.The specific impulse is a measure of how effectively a propulsion system uses a propellant.By definition, it is total impulse delivered per unit of propellant.The impulse bit is the impulse value of every pulse time, and efficiency is defined as the ratio of kinetic energy to discharge energy.Solving Equation ( 29), we can obtain the velocity and displacement of the current sheet, the discharge voltage, and the current.On this basis, these performance parameters can be calculated further as follows [14,19]: .
x(t ) m bit g = .
x(t ) g where I sp stands for the specific impulse, I bit is the impulse bit, η stands for efficiency, T/P is the ratio of thrust to power, t is the time of a single pulse, and E is the consumed energy at every pulse time.

Method Validation
In order to validate the new modeling method, we tested the accuracy of the model and the assumptions throughout by comparing the simulation discharge waveforms, heat flux, and mass flux with experimental data from other papers [12,19,31,32].Moreover, for further comparison with the simulation results, we adopted some experiments to measure performance parameters except specific impulse, impulse, and single pulse ablative mass.
Parameters of Teflon under different states are shown in Table 1.

Discharge Waveform Comparison
The simulation discharge waveforms were compared with data from two parallel plate thrusters: the µPPT in Coletti's paper [19] and the LES-6 PPT in Vondra's papers [12,32].Table 2 gives the operation parameters of µPPT by Coletti [19], and Table 3 gives the operation parameters of LES-6 PPT by Vondra [12].Figure 3 presents comparisons of discharge current and voltage waveforms at different discharge energies between the simulation results and reference data from Coletti [19].Similar comparison with the reference data from Vondra [12,32] are shown in Figure 4.By comparing these graphs, we see that the simulation results are in agreement with the experimental results; some disagreements are mainly attributed to the assumption that the thruster is under an ideal and stable working state.In the simulation process, the thruster is assumed to have always been working steadily.However, in the actual working process of the thruster, it needs a lag time to enter the stable working state.This unstable operation results in disagreement at the start of the discharge process, and the lag time leads to a time shift between the simulated and reference data.Other mismatch mainly results from the assumption that the plasma resistance is constant.For the simulation process, the value of the plasma resistance was calculated by Equation ( 28) to be a constant value.But in the real thruster, the plasma resistance varies as the plasma is being formed and as the electron temperature and number density change [15].Therefore, the simulation results agreed better with the experimental results, as the plasma resistance had not changed much in the early discharge process.With the discharge process going on, the variation of plasma resistance caused some disagreement between the simulation results and the experimental results.

Ablation Parameters Comparison
To further validate the reliability of the model, especially the ablation process, we compared the simulation results of the heat flux and mass flux with the data from Stechman's thesis [31].The numerical results agree with the given experimental data in permissible error range, as presented in Figure 5.The disagreement between the experiment and the simulation may be attributed to the absorption ratio of the discharge energy (γ) selected in the simulation process.In the simulation process, the ratio of single pulse discharge energy absorbed by the surface of the Teflon solid (γ) was assumed to be 0.05.But in the real thruster, this value will fluctuate slightly at different times rather than always being constant.

Perfomance Parameters Comparison
This paper adopts the torsional pendulum built in [33] to measure the impulse bit.A photograph and schematic illustration of the thrust stand are shown in Figure 6.As shown in Figure 6a, the torsional pendulum inside the vacuum tank consists of the following: (1) a pulsed plasma thruster, (2) a torsional arm, (3) a Csection torsional bar made of beryllium-bronze, (4) a circular mirror, (5) a He-Ne laser, (6) rectangular mirrors (with silver plating), (7) a counterweight, (8) a one-dimensional position sensing detector (PSD) to measure the micro displacement of the He-Ne laser spot, (9) an electromagnetic calibration device consisting of a pulsed ampere force generator inside the vacuum tank and a function generator outside the vacuum tank, and (10) an electrical adapter for the electrical connection between the thruster inside the vacuum tank and the power supply and ignition controller outside the vacuum tank.As shown in Figure 6b, the pulsed plasma thruster is set onto the torsional pendulum.When the thruster works, the impulse it produces makes the torsional arm have a jigging motion.Additionally, the swing of the torsional arm drives the circular mirror fixed on the arm to swing together.The laser beam irradiates the circular mirror and eventually the PSD sensor, after being reflected several times by mirrors.Following this, the location of the laser beam spot is recorded in real time by the PSD.Fixing the laser and the reflected laser beam is deflected correspondingly when the circular mirror swings.Under the range of the torsional pendulum, the maximum deviation from the balance position of the laser spot on the PSD is proportional to the impulse value of the arm.Thus, the impulse of the thruster can be calculated as well [33].The electronic balance was used to calculate the ablation mass.Measuring the mass variation value of the Teflon before and after the multiple pulse discharge process, the average ablation mass in single pulse time can be obtained.
The selected operation parameters in our experiment were the same as those in Table 2. Following this, the performance parameters obtained in the experiment are shown together with the simulation results in Table 4.The experimental results in Coletti's paper [19] are also attached to Table 4 to validate the simulation results and experimental results obtained by our experiments.The numerical results agree with the given experimental data in the permissible error range, as presented in Table 4. Furthermore, we also used the operation parameters in Table 3 for the simulation performance parameters, which were directly compared with the experimental data of the LES-6 PPT in Zhang's paper [14].The simulation results are in agreement with the experimental results, as given in Table 5.The results discussed above prove the reliability of the model, which is also less computation-intensive.Not only was rapid prediction of the µPPT macroscopic performance realized, but the influences of different variables on µPPT performance could be explored with this simulation method, which also provided a guide for the design of the high-performance µPPT.

Parameter Analysis
Parameters that affect µPPT performance mainly include electrical parameters, structural parameters, and propellant types.The electrical parameters mainly include circuit capacitance, inductance, resistance, charging voltage, and discharge energy, and the structural parameters mainly take into account plate configuration, working mass supply mode, and spark plug structure and arrangement.In this section, if there is no special description, the operation parameters adopted are the same as the ones presented in Table 2.
The use of an experimental method to explore the effects of all these factors is time-consuming and wastes money.With this validated modeling method, these parameters can be analyzed accurately in a short time.

Initial Discharge Voltage
The single pulse discharge energy increases with variation of the initial discharge voltage of the capacitor when the other parameters remain unchanged.Figure 7 gives the simulation results at different initial discharge voltages.
The peak value of the discharge current and voltage increases significantly with increased initial discharge voltage, while the discharge cycle, which refers to the time interval between two adjacent peaks or two adjacent troughs in the discharge waveforms, remains constant, as shown in Figure 7a,b.It is worth noting that oscillation would occur with increased initial discharge voltage, which exacerbates the impact on the capacitor and causes the capacitor to lose efficacy.
Figure 7c shows variation of the surface ablation temperature: when the initial voltage increases, the surface ablation temperature goes up, so the second ablation stage lasts longer.As a result, the rate of ablation is accelerated and the corresponding ablation mass increases, as shown in Figure 7d.
Figure 7e shows that when the initial value of the discharge voltage increases, the impulse bit becomes larger, while the specific impulse decreases gradually, which is mainly due to the increased single pulse ablation mass.Efficiency has a positive correlation with initial discharge voltage, while the thrust-power ratio is just the opposite, as shown in Figure 7f.

Capacitance
The single pulse discharge energy also has a positive correlation with capacitance.The capacitance values were 1 µF, 2 µF, and 5 µF, and the numerical results are given in Figure 8.
As shown in Figures 6b and 8a, with increased capacitance, the peak of the discharge current and voltage becomes larger, while the discharge time becomes longer.However, the reverse voltage and current decrease, which reduces the impact of the current on the capacitor.
Larger capacitance means larger discharge energy and smaller fluctuation of the discharge waveform.Variation of the surface ablation temperature is shown in Figure 8c; it rises with the increase of the capacitor, while the duration when the temperature is higher than the melting temperature is longer.Additionally, following this, the ablation mass accumulates.
Figure 8e,f show that the impulse bit, the specific impulse, and the efficiency increase with increased capacitance.However, the thrust power ratio has a downward trend with increased capacitance.The main reason for this may be that the power increases as the capacitance increases, and that the improvement of capacitance on the power is greater than that on the impulse bit; this agrees with the investigation of Schonherr [34].
In summary, larger capacitances can effectively reduce the circuit oscillation and improve the comprehensive performance.However, to raise the capacitance, its whole mass and volume will become larger, which needs to be taken into consideration with regard to the overall mission requirements.In order to explore the effects of different capacitance-initial voltage ratios on the comprehensive performance of µPPTs, the single pulse discharge energy was set at 3.84 J and the corresponding capacitance and voltage values were 1 µF/2771 V, 3 µF/1600 V, and 5 µF/1239 V.The simulation results are given in Figure 9.
The peak value of the discharge current and ablation surface temperature decreases, but the specific impulse and efficiency improve when the initial discharge voltage is reduced and the capacity of the capacitor increases under the same discharge energy.What is more, the degree of oscillation of the discharge wave is weakened, which helps to reduce the impact of the current on the capacitor.As a result, the ablation mass is reduced while the impulse bit and the thrust-power ratio are basically invariant.

External Resistance
External resistance, which usually includes capacitance resistance, transmission line resistance, and plate resistance, is one of the important factors that cause µPPT discharge energy loss.For the simulation calculation, the external resistances were selected to be 3 Ω, 33 Ω, and 63 Ω, keeping the other parameters constant.The results are presented in Figure 10.The values of the resistances were chosen randomly.
As we can see in Figure 8, the crest value of the voltage and current decreases with increased resistance.Meanwhile, the discharge cycle remains roughly constant and the oscillation becomes smaller.As far as ablation surface temperature is concerned, its peak value is insensitive to resistance and time, while the time during which the temperature is above melting temperature becomes shorter.The ablation mass, specific impulse, impulse bit, efficiency, and thrust-power ratio decrease, mainly due to the Ohmic heat loss and the shorter ablative process.In short, increasing resistance will exacerbate energy loss and seriously affect the comprehensive performance; large external resistance should therefore be avoided.

Plate Gap
The plate gap directly affects the distribution of plasma resistance and arc on the ablation surface, which will also have an important influence on the comprehensive performance and ignition reliability.
The plasma resistance and the ablation surface will also become larger with increased plate gap.Regarding these factors, the extreme values of the waveforms and the degree of fluctuation decrease, as shown in Figure 11a,b, while the surface temperature decreases and the time during which the temperature is higher than the melting temperature also decreases, as shown in Figure 11c.Also, the ablation mass decreases due to the shorter second ablation stage, as seen in Figure 11d.
What is more, the specific impulse, impulse bit, thrust-power ratio, and efficiency are improved, but there are still some differences.If the plate is in a relatively large space, it can improve the specific impulse and efficiency by increasing the distance more validly.Meanwhile, if the plate is in a relatively small space, the impulse bit and thrust-power ratio can be improved more validly, as expressed in Figure 11e,f.

Plate Width
Enlarging the plate width reduces the plasma resistance while increasing the ablation surface and then lowering the surface energy density.As can be seen in Figure 12, when the plate width increases, the crest value of the discharge waveforms rises, the fluctuation degree becomes more intense, and the peak value of the discharge current increases, but the surface temperature declines a bit.The ablation mass, impulse bit, and efficiency are in negative relation to the plate width, while the specific impulse exhibits an opposite trend.

Different Plate Gap-Width Ratios under the Constant Ablation Area
In order to explore the influence of different gap-width ratios on the working process of µPPTs, the single pulse discharge energy was set at 3.84 J, and the corresponding values of plate gap and width were 10 mm/20 mm, 14.14 mm/14.14mm, 20 mm/10 mm, and 24.5 mm/8.16mm.
The peak value of the discharge current increases significantly, while the crest value of discharge voltage declines with increased initial discharge voltage and the cycle remains steady, as shown in Figure 13a,b.Figure 13c shows that the maximum ablation temperature keeps still, but the time during which the temperature is higher than the melting temperature decreases.The ablation mass is negatively correlated with the plate gap-width ratio, as seen in Figure 13d.
It can be seen in Figure 13e,f that the comprehensive performance became better with an increased ratio.

Conclusions
This paper introduces a modeling method that combines an ablation model with a modified electromechanical µPPT model for a more accurate µPPT simulation.The model adopts an integration method to reflect the mass variation in the discharge process, which differs from other modeling methods.What is more, the non-Fourier effect during the ablation process was also taken into consideration.
The numerical results are validated against the experimental values, and good agreement is found.Furthermore, the effects of different operation parameters are explored using the developed model, and the main conclusion can be drawn as follows: 1.
Different ways of increasing the discharge energy have different effects on the performance of µPPT.Larger capacitances can effectively weaken the oscillation characteristics of the circuit and improve the performance, except the impulse bit, specific impulse, and efficiency.Meanwhile, an increase in the initial discharge voltage enlarges the impulse bit but reduces the specific impulse due to the explosion of the ablation mass in the high initial discharge voltage.Moreover, under constant discharge energy, the combination of high capacitance and small initial voltage can enhance the specific impulse.

2.
The loss of Ohmic heat, which is due to large resistance, causes the loss of discharge energy, resulting in decreased performance, including the specific impulse, impulse bit, and efficiency.Therefore, larger circuit resistances should be avoided in order to prevent the decrease in performance from becoming serious.

3.
Under the constant ablation surface area, the ablation mass, peak value of discharge current, and degree of oscillation decline with increased plate gap-width ratio.On the contrary, the impulse bit, specific impulse, and efficiency have a positive relation with the plate gap-width ratio.Thus, a high plate gap-to-width ratio can improve µPPT performance.

Figure 2 .
Figure 2. Schematic of the physical model and coordinates.(a) The first ablation stage; (b) The second ablation stage; (c) The third ablation stage.

Figure 3 .
Figure 3. Comparisons of discharge current and voltage waveforms at different discharge energies between simulation results and reference data from Coletti [19]: (a) at 5 J discharge energy; (b) at 4 J discharge energy; and (c) at 3 J discharge energy.

Figure 4 .
Figure 4. Comparison of discharge current and voltage waveforms between simulation results and reference data from Vondra [12,32].

Figure 5 .
Figure 5.Comparison of heat flux and mass flux between simulation results and reference data from Stechman [31]: (a) heat flux variation with time; and (b) mass flux variation with time.

Figure 7 .
Figure 7. Simulation results at different initial discharge voltages: (a) Discharge voltage waveform; (b) discharge current waveform; (c) variation of ablation surface temperature with time; (d) variation of ablation mass with initial discharge voltage; (e) variation of impulse bit and specific impulse with initial discharge voltage; and (f) variation of efficiency and ratio of thrust to power with initial discharge voltage.

Figure 8 .
Figure 8. Simulation results at different capacitances: (a) Discharge voltage waveform; (b) discharge current waveform; (c) variation of ablation surface temperature with time; (d) variation of ablation mass with capacitance; (e) variation of impulse bit and specific impulse with capacitance; and (f) variation of efficiency and ratio of thrust to power with capacitance.

Figure 9 .
Figure 9. Simulation results at different capacitance-initial voltage ratios: (a) Discharge voltage waveform; (b) discharge current waveform; (c) variation of ablation surface temperature of with time; (d) variation of ablation mass with capacitance; (e) variation of impulse bit and specific impulse with capacitance; and (f) variation of efficiency and ratio of thrust to power with capacitance.

Figure 10 .
Figure 10.Simulation results at different external resistances: (a) Discharge voltage waveform; (b) discharge current waveform; (c) variation of ablation surface temperature with time; (d) variation of ablation mass with external resistance; (e) variation of impulse bit and specific impulse with external resistance; (f) variation of efficiency and ratio of thrust to power with external resistance.

Figure 11 .
Figure 11.Simulation results at different plate gaps: (a) Discharge voltage waveform; (b) discharge current waveform; (c) variation of ablation surface temperature with time; (d) variation of ablation mass with plate gap; (e) variation of impulse bit and specific impulse with plate gap; (f) variation of efficiency and ratio of thrust to power with plate gap.

Figure 12 .
Figure 12.Simulation results at different plate widths: (a) Discharge voltage waveform; (b) discharge current waveform; (c) variation of ablation surface temperature with time; (d) variation of ablation mass with plate width; (e) variation of impulse bit and specific impulse with plate width; (f) variation of efficiency and ratio of thrust to power with plate width.

Figure 13 .
Figure 13.Simulation results at different plate gap-width ratios: (a) Discharge voltage waveform; (b) discharge current waveform; (c) variation of ablation surface temperature with time; (d) variation of ablation mass with ratio; (e) variation of impulse bit and specific impulse with ratio; (f) variation of efficiency and ratio of thrust to power with ratio.

Table 4 .
Comparison of performance parameters between experiment and simulation.

Table 5 .
[14]arison of performance parameters between simulation results and experimental data of the LES-6 PPT in Zhang's paper[14].
Author Contributions: Y.O.conceived and designed the experiments; J.W. and S.T. performed the experiments; Y.O. and Y.Z.analyzed the data; J.L. contributed reagents/materials/analysis tools; Y.O.wrote the paper.The authors declare no conflict of interest.Q m energy released per unit time and volume, W/m 3 R universal gas constant, Ω R c resistance of the capacitor, Ω R e resistance of the wires and leads, Ω R ps resistance of the plate electrodes, Ω Acknowledgments: This work was supported by Chinese National Natural Science Foundation (11772354).Conflicts of Interest: