2.1. Research objective

In this paper, the research objective is the ECRIT developed by the Harbin Institute of Technology faced with the space GW detection task. The structure and a photo of the thruster are illustrated as figure 1. The propellant Xe is fed into the discharge chamber through a flow tube and a gas distributor. The MW is injected into the thruster through an antenna. The opposite direct current (DC) voltages are applied to the SG and AG, which can form a strong electric field between the grids and accelerate ions to produce thrust. Besides, it adopts two-stage permanent magnets to form the cusped magnetic field and confine the electrons. On one hand, the magnetic configuration can enhance ionization intensity and improve specific impulses. On the other hand, the magnetic configuration can reduce the interaction between plasma and walls, which can increase thruster lifetime. The diameter of the discharge chamber is 1 cm. There are 85 holes of grids. The diameters of SG holes and AG holes are respectively 0.72 mm and 0.4 mm. The thicknesses of SG and AG are respectively 0.2 mm and 0.3 mm. The distance between two grids is 0.25 mm. Experimental results show its thrust range is 1∼100 μN by adjusting operating parameters. The ERCIT should be ignited before the detection task starts. Under the detection mode, the drag-free controller generates thrust commands according to the non-conservative forces on satellites. The operating parameters of ECRITs are automatically and continuously adjusted by the feedback controller.

Figure 1. The illustration of ECRIT. (a) Structure and (b) photo.

2.2. Experimental set-up

Experiments are conducted to obtain necessary data for the simulation model and its validation. The experimental platform is illustrated as figure 2. The vacuum chamber is cuboid, whose size is 0.8 × 0.8 × 0.7 m. Two molecular pumps and one mechanical pump work to maintain the pressure below 10⁻⁴ Pa. The mass rate of propellants is controlled by an Alicat flow meter. The DC sources are developed by HSPY, whose ranges are 0∼1500 V. The MW power sources can realize a 0∼20 W range. These devices satisfy the experimental requirements.

Figure 2. Experimental platform.

There are two experiments which are conducted in the paper. Firstly, the electron energy distribution needs to be tested by a Langmuir probe because it deviates from a Maxwell distribution as it is affected by the complex cusped magnetic field and many inelastic collision processes. The real electron energy distribution is critical for computing the rate coefficients of the ionization and excitation reactions. The Langmuir probe consists of a cylindrical tantalum and a ceramic tube, whose diameter is 0.7 mm. Considering the strongly magnetized plasma in the discharge chamber, the planar probe is utilized to avoid some problems of cylindrical probes like radical sheath expansion and electron non-saturation. The flat end of the Langmuir probe is located at the ECR layer, which is 4.5 mm from the SG axially and 2 mm from the central axis radically. The scanning range of probe voltage is −35–150 V. The electron energy distribution is obtained by differentiating the transition region of the I–V characteristic curve, which is shown in § 2.3. Secondly, we measure SG current and AG current under different operations to validate the model accuracy. These electric parameters are recorded by a YOKOGAWA wave recorder. These experimental results are shown at the beginning of §§ 3.1–3.3.

2.3. Simulation model

The integrative model is based on some assumptions listed as follows. The plasma parameters in the discharge chamber are uniform; the sheath profile upstream of the SG is approximately elliptic paraboloid (Yan et al. Reference Yan, Guo, Wu, Gu and Yang2020); the ion beamlet shape is approximately a rotating hyperboloid considering the Coulomb repulsion effect of ions in the grid system (Shang et al. Reference Shang, Xiang, Jiang, Wang, He and Wen2019; Li et al. Reference Li, Yang, Zhang, Zhang and Sun2022). The main structure of the model is shown in figure 3. The plasma parameters in the discharge chamber are simulated by a global model. The discharge chamber model is connected with the sheath model by plasma parameters upstream of the SG and grid voltage. The ion beamlet model is connected with the sheath model by sheath profile parameters.

Figure 3. The structure of the model.

It can be seen from figure 3 that the influences of the flow rate and MW power on the ion acceleration characteristics are considered via the sheath potential in this model. Besides, the influence of grid voltage on the ion extraction properties is considered via variations of the sheath shape in the model. Compared with traditional global model, the integrative model built in this paper can not only analyse plasma parameters in the discharge chamber, but also analyse variations of the sheath shape with operating parameters. Besides, this model considers the divergence of the ion beamlet and the intersection effect of the AG, which are helpful for analysing thruster dynamic processes with adjustment of operating parameters. The main equations and parameters can be referenced in Niu, Liu & Yu (Reference Niu, Liu and Yu2024), which are introduced briefly as follows. The discharge chamber model includes the ion number conservation equation (Lafleur et al. Reference Lafleur, Habl, Rossi and Rafalskyi2022), the atom number conservation equation (Magaldi et al. Reference Magaldi, Karnopp, Sobrinho and Pessoa2022), the electron energy conservation equation (Lucken et al. Reference Lucken, Marmuse, Bourdon, Chabert and Tavant2019; Nakagawa et al. Reference Nakagawa, Koizumi, Naito and Komurasaki2020) and the atom energy conservation equation (Chabert et al. Reference Chabert, Arancibia Monreal, Bredin, Popelier and Aanesland2012; Grodein et al. Reference Grodein, Lafleur, Chabert and Aanesland2016). These equations are solved to update the ion density, electron temperature, atom density and atom temperature.

The rate coefficients of the ionization reaction and excitation reaction are two key parameters to calculate the plasma parameters in the discharge chamber, which are determined by the cross-sections of collisions and the electron energy distribution function. The cross-sections of the ionization and excitation reaction are as in Niu et al. (Reference Niu, Liu and Yu2024). However, the electron energy distribution function deviates from a Maxwell distribution affected by the complex cusped magnetic field and many inelastic collision processes (especially the ionization and excitation process). Therefore, the ECRIT is diagnosed by a Langmuir probe in the ECR region to obtain its electron energy distribution function.

The diagnostic result of the electron energy distribution is obtained by differentiating the transition region of I–V characteristic curves. Several typical electron energy distribution functions with parameters, shown as (2.1), are combined to approach the measurement results. Here, ${T_e}$ is equivalent electron temperature. The computations of ${c_{1i}}$ and ${c_{2i}}$ are shown as (2.2). Parameter i is the characteristic parameter of the distribution function. We can obtain a Maxwell distribution when i equals to 1 and a Druyvesteyn distribution when i equals to 2. Here, $\varGamma$ is gamma function. The parameters are fitted by the least square method and then revised to satisfy the definition, shown as (2.3). The diagnostic and fitted results are shown in figure 4(a). The rate coefficients of ionization and excitation reactions are calculated by (2.4) and shown in figure 4(b). Here, Q and $\sigma (E)$ are respectively the rate coefficient and cross-section of corresponding collision process and m is the electron mass

(2.1)\begin{gather}f(E) = \sum\limits_{i = 1}^n {{c_{1i}}T_e^{ - 1.5}\sqrt E \exp ( - {c_{2i}}{{(E/{T_e})}^i})} ,\end{gather}

(2.2)\begin{gather}\left. {\begin{array}{c@{}} {{c_{1i}} = i{{\left( {\dfrac{2}{3}} \right)}^{3/2}}\dfrac{{{{[\varGamma (5/2i)]}^{3/2}}}}{{{{[\varGamma (3/2i)]}^{5/2}}}}}\\ {{c_{2i}} = {{\left( {\dfrac{2}{3}} \right)}^{3/2}}{{\left[ {\dfrac{{\varGamma (5/2i)}}{{\varGamma (3/2i)}}} \right]}^i}} \end{array}} \right\},\end{gather}

(2.3)\begin{gather}\int_0^\infty {f(E)\,\textrm{d}E} = 1,\end{gather}

(2.4)\begin{gather}Q = \int_0^\infty {\sigma (E)} \sqrt {\frac{{2E}}{m}} f(E)\,\textrm{d}E.\end{gather}

Figure 4. Electron energy distribution function and the calculated rate coefficient. (a) Electron energy distribution function and (b) reaction rate coefficient.

The absorbed MW power is shown as (2.5) and (2.6) (Budden Reference Budden1966). To reflect the saturation effect of MW power absorption, the model sets the upper limit of absorbed MW power when the plasma density reaches its critical value, shown as (2.7). The relationship between the SG sheath potential and electron temperature upstream of the SG is shown as (2.8) (Goebel & Katz Reference Goebel and Katz2008)

(2.5)\begin{gather}{P_{\textrm{abs}}} = {P_{\textrm{input}}}{\eta _p},\end{gather}

(2.6)\begin{gather}{\eta _p} = (1 - {\textrm{e}^{ - \mathrm{\pi }\eta }}),\end{gather}

(2.7)\begin{gather}\eta = \frac{{\omega _{\textrm{pe}}^2}}{{\omega c\alpha }};\quad {\omega _{\textrm{pe}}} = \sqrt {\frac{{{n_{\textrm{ecr}}}{e^2}}}{{\varepsilon m}}} ;\quad \alpha = \frac{1}{B}\frac{{\textrm{d}B}}{{\textrm{d}z}},\end{gather}

(2.8)\begin{gather}{n_c} = {({f_p}/9)^2},\end{gather}

(2.9)\begin{gather}{\phi _e} = {T_e}\left\{ {\ln \left( {\sqrt {\frac{M}{{2\mathrm{\pi }m}}} } \right) - \ln \left( {\frac{{{A_i}}}{{{A_e}}}} \right)} \right\}.\end{gather}

Here,${P_{\textrm{input}}}$ is the input MW power; ${\eta _p}$ is the absorbed efficiency of plasma for MW power; ${f_p}$ is the MW frequency, 4.2 GHz; c is the light speed; $\alpha$ is the normalized magnetic induction intensity gradient in the ECR region; ${\omega _{\textrm{pe}}}$ is the plasma oscillation angular frequency; ${n_{\textrm{ecr}}}$ is the average ion density in the ECR region, satisfying ${n_{\textrm{ecr}}} = {f_{\textrm{ecr}}}n$, n is the average electron density in the discharge chamber; $\varepsilon$ is the permittivity of free space; ${n_c}$ is the critical plasma density considering the saturation effect of MW power absorption by plasma. Also, M and m are respectively the Xe atom and electron masses and ${A_i}$ and ${A_e}$ are respectively the ion and electron loss areas on the SG.

We adopt the elliptic paraboloid to describe the sheath profile upstream of the SG considering the radial and angular symmetry of ECRITs, shown as (2.9). The origin of the coordinate system is located at the intersection point of the sheath and the central axis of the thruster. The z-axis is the direction of the central axis of the thruster. The x-axis and y-axis are two orthogonal axes perpendicular to the z-axis, satisfying the right-hand spiral rule. The ion beamlet shape is described by a rotating hyperboloid, shown as (2.10). Symbol ‘b ₁’ is the minimum cross-sectional area of the ion beamlet and ‘a’ is its z-direction location; symbol ‘d ₀’ describes the bending degree of the ion beamlet. The relationships between these parameters with sheath profiles are shown in Niu et al. (Reference Niu, Liu and Yu2024)

(2.10)\begin{gather}{x^2} + {y^2} = 2pz,\end{gather}

(2.11)\begin{gather}{x^2} + {y^2} = {d_0}{(z - a)^2} + {b_1}.\end{gather}

The reference value of the sheath thickness is the distance from the front side of the SG to the centre of the AG, shown as (2.11). Symbols ${l_s}$ and ${l_a}$ are respectively the thicknesses of the SG and the AG. Symbol ${l_{sa}}$ is the distance between the SG and the AG. The sheath thickness L_s is the distance from the origin to the centre of the AG, shown as (2.12) and related to the grid voltages, sheath potential and plasma density upstream of the SG. Symbols ${f_{\textrm{ls}}}$ and $\zeta$ are parameters of the sheath thickness. Symbol ${\lambda _d}$ is the Debye length of plasma upstream of the SG. The sheath profile is determined by the relationship between L_s and L_{s 0}, which is parallel to the SG when L_s equals to L_{s 0}

(2.12)\begin{gather}{L_{s0}} = {l_s} + {l_{sa}} + 0.5{l_a},\end{gather}

(2.13)\begin{gather}{L_s} = {f_{\textrm{ls}}}{\lambda _d}{\left( {\frac{{e({\phi_{\textrm{es}}} + {U_s} - {U_a})}}{{k{T_{\textrm{es}}}}}} \right)^\zeta }.\end{gather}

We choose a single-aperture model to study the sheath and ion beamlet characteristics with the assumption that the plasma density is uniform upstream of the SG. The simulation object is a cylinder, whose radius and height are respectively the radius of the SG aperture and the distance from the front sidewall of the SG to the back sidewall of the AG. The front sidewall refers to the surface close to the discharge chamber. The ion number conservation equation in the grids is shown as (2.13). The term on the left side of (2.13) is the temporal gradient of the average ion density. The terms on the right side of (2.13) are respectively the ion density change rates caused by ions entering the SG, extracted from the AG and impinging on the AG (sidewall and inner wall of the AG aperture). The calculation of parameters in (2.13) are shown in Niu et al. (Reference Niu, Liu and Yu2024)

(2.14)\begin{equation}\frac{{\textrm{d}{n_v}}}{{\textrm{d}t}} = \frac{{{n_s}{A_s}{T_{s2}}\sqrt {\frac{{e{T_{\textrm{es}}}}}{M}} }}{{N{V_{\textrm{grid}}}}} - \frac{{{f_{v1}}{n_v}{\eta _{f1}}{A_{\textrm{ah1}}}{v_{i1}}}}{{{V_{\textrm{grid}}}}} - \frac{{({f_{v2}}{\eta _{f2}}{A_{\textrm{ah2}}} + {f_{v3}}{\eta _{f3}}{A_{\textrm{ah3}}}){n_v}{v_{i2}}}}{{{V_{\textrm{grid}}}}}.\end{equation}

Here, ${n_v}$ is the average ion density for the single-aperture body; N is the SG aperture number; ${V_{\textrm{grid}}}$ is the volume of single-aperture body; ${f_{v1}}$, ${f_{v2}}$ and ${f_{v3}}$ are respectively the ratio of ion density on the AG exit, AG sidewall, AG aperture inner wall and ${n_v}$; ${\eta _{f1}}$, ${\eta _{f2}}$ and ${\eta _{f3}}$ are the area integrals of the cosine of the angle respectively between the extracted ion acceleration direction, the ion (impinging on the AG sidewall) acceleration direction, the ion (impinging on the AG aperture inner wall) acceleration direction and the z direction; ${A_{\textrm{ah1}}}$ is the extraction area of the ion beamlet; ${A_{\textrm{ah2}}}$ is the impingement area on the AG sidewall; ${A_{\textrm{ah3}}}$ is the impingement area on the AG aperture inner wall; ${v_{i1}}$ and ${v_{i2}}$ are the speeds of ions respectively extracted from the AG and impinging on the AG.

The SG current, AG current, AG current from charge exchange (CEX) ions and thrust are illustrated as (2.14)∼(2.17). The AG current includes three parts, namely the primary ions impinging on the AG sidewall and AG aperture inner wall from the discharge chamber and the CEX ions impinging on the AG sidewall from downstream of the AG (Lingwei et al. Reference Lingwei, Yu, Juan, Zuo, Haocheng, Haixing and Haibin2010). The CEX part is shown as (2.16), whose parameters can be referenced in Niu et al. (Reference Niu, Liu and Yu2024). The thrust includes two parts, which are respectively produced by ions and neutral gas. The atom speed is shown as (2.18). Symbols ${v_a}$ and ${T_g}$ are connected by the ideal gas state equation. The unit of ${T_g}$ is eV, which is updated by the atom energy conservation equation

(2.15)\begin{gather}{I_b} = {n_s}e{A_s}{T_{s2}}\sqrt {\frac{{k{T_{\textrm{es}}}}}{M}} ,\end{gather}

(2.16)\begin{gather}{I_a} = {f_{v2}}{n_v}e{\eta _{f2}}{N_s}{A_{\textrm{ah2}}}{v_{i2}} + {f_{v3}}{n_v}e{\eta _{f3}}{N_s}{A_{\textrm{ah3}}}{v_{i2}} + {I_{\textrm{cex}}},\end{gather}

(2.17)\begin{gather}{I_{\textrm{cex}}} = \int_0^{{z_0}} {{f_I}e{A_a}\frac{{\textrm{d}{n_{\textrm{cex}}}}}{{\textrm{d}t}}\,\textrm{d}z} ,\end{gather}

(2.18)\begin{gather}F = {\eta _{f1}}\frac{{({I_b} - {I_a})}}{e}\sqrt {\frac{{2M({\phi _{\textrm{es}}} + {U_s})}}{e}} + \frac{1}{4}{n_g}\sqrt {\frac{{8e{T_g}}}{{\mathrm{\pi }M}}} {A_a}{T_a}{\eta _c}M{v_a},\end{gather}

(2.19)\begin{gather}{v_a} = \sqrt {3e{T_g}/M} .\end{gather}

Here, ${I_b}$ is the SG current; ${n_s}$ is the plasma density upstream of the SG; ${A_s}$ is the SG area; ${T_{s2}}$ is the effective transparency of the SG; ${T_{\textrm{es}}}$ is the electron temperature upstream of the SG; ${I_a}$ is the AG current; ${\phi _{\textrm{es}}}$ is the sheath potential upstream of the SG; ${U_s}$ is the SG voltage; ${T_a}$ is the optical transparency of the AG; ${\eta _c}$ is the Clausing factor of the grid system; ${I_{\textrm{cex}}}$ is the AG current due to CEX reaction downstream of the AG; ${n_{\textrm{cex}}}$ is the density of CEX ions; ${f_I}$ is the factor describing the proportion of CEX ions hitting the AG; ${v_a}$ is the atom speed; ${T_g}$ is the neutral gas temperature.

In this paper, the integrative model is built to study the influences of operating parameters on the dynamic process of ECRITs. Therefore, the model accuracy needs to be validated. Firstly, the stable accuracies of the model with variations of MW power, flow rate and grid voltages are all validated in § 3 by experiments. Besides, the dynamic accuracy of the model is validated by adjusting MW power. This is mainly because the response time of the MW source is approximately 200 μs, much shorter than that of flow rate controller and DC sources used in the experiments. Therefore, the dynamic characteristics of ECRITs can be reflected from the grid current by adjusting MW power. In § 3, there is an analysis about the response time of the plasma parameters which refers to their stable time after operations are adjusted.

3.1. Influence of MW power on thruster dynamic characteristics

When the flow rate and SG/AG voltages are set as 0.03 sccm and 600 V/−150 V, the SG current under different MW powers is shown in figure 5(a). Sccm is the unit of flow rate representing standard millilitres per minute. It can be seen that the simulation results agree well with experiments, with an error less than 5 %. Then the MW power is adjusted from 0.15 to 0.6 W by the MW source control system while keeping the flow rate and grid voltages invariant. The MW power is measured by the wave detector and its response time is approximately 200 μs. We choose this process as the model input and obtain the corresponding simulation results, which agree well with experiments, shown in figure 5(b). It shows that the response time is approximately 10 ms, which is much longer than the response process of the MW power and can reflect the dynamic characteristics of ECRITs with the adjustment of MW power.

Figure 5. Comparison of experiments and model results with the variation of MW power. (a) Stable characteristics and (b) dynamic characteristics.

Experimental results show that the characteristic time of the ECRIT dynamic process is of the order of 10 ms if the MW power can be adjusted much more rapidly. Actually, the characteristic time of the MW transfer velocity and electron acceleration by cyclotron resonance determines that the energy can be absorbed by electrons instantaneously after the adjustment of the MW power. However, rapid increase of the electron temperature breaks the ionization balance, which results in the motion of neutral gas and redistribution of atom density. Considering the volume of the discharge chamber and atom motion velocity, the recovery time of ionization balance is approximately 10 ms. More detailed analysis is listed as follows.

The flow rate and SG/AG voltages are respectively set as 0.03 sccm and 600 V/−150 V. The state parameters in the discharge chamber are shown as figure 6 when the MW power steps from 0.15 to 0.45 W with 0.15 W intervals. In the figures showing the dynamic processes of thruster parameters with the adjustment of operating parameters, we use the blue dashed line to represent the adjustment of operating parameters (MW power, flow rate and DC voltage) and use the red solid line to represent the dynamic processes of thruster parameters. It can be seen from figure 6(a) that the electron temperature firstly increases with enhancement of the MW electric field after the MW power rises. The increase of high-energy electrons can enhance the ionization process and improve electron density. However, the ionization process is a strong energy dissipation process. Therefore, the promotion of ionization intensity costs more electron energy and the electron temperature falls with the increase of electron density. The response time of this stage is tens of microseconds, which is consistent with the plasma characteristic time.

Figure 6. Simulation results of state parameters with step response of MW. (a) Electron temperature, (b) electron density, (c) atom density, (d) atom temperature, (e) source term of collisions of ions and atoms and (f) source term of collisions of ions and wall.

The response time of the next stage is approximately 10 ms due to the volume effect. The enhancement of the ionization process causes a drop of the atom density locally, especially in the ECR region, and further results in the diffusion motion of neutral gas. The motion of neutral gas causes the decline of average atom density of the discharge chamber, including the main ionization region (the cusp region). The decrease of the average atom density further results in the fall of ionization intensity, which can be reflected from the electron density at this stage. Therefore, the average electron density firstly increases and then decreases, which forms an overshoot with the step increase of MW power. The actual volume of the discharge chamber determines the time scale of neutral gas filling and it is of the order of tens of milliseconds. Therefore, the volume effect can be reflected from the atom density in figure 6(c).

It can be seen from figure 6(d) that the response process of atom temperature can also be divided into two stages. The response time of the first stage is approximately 0.2 ms, longer than the plasma response time with the increase of MW power. The main processes influencing atom temperature include the collision process between ions and neutral gas and the collision process between neutral gas and the discharge chamber wall. The source terms of these two processes are shown as figures 6(e) and 6(f). It can be seen that the characteristic time of the atom temperature variation due to collisions between ions and neutral gas is tens of microseconds, consistent with the plasma response. The characteristic time of the atom temperature variation due to collisions between neutral gas and the wall is 0.2 ms, which is the main factor determining the response time of atom temperature in the first stage.

The SG current and thrust with a step increase of MW power are shown as figure 7. Under the adjustment process, the thruster is optimal perveance and the AG current is very low. Therefore, the thrust is mainly determined by the SG current with grid voltages unchanged. The SG current firstly increases after MW power increases, which characteristic time is tens of microseconds caused by the enhancement of the ionization process. The characteristic time of the next stage is approximately 10 ms, which is the characteristic time of neutral gas diffusion motion resulting from the broken ionization balance. The overshoots of the SG current and thrust occur because of the decrease of average atom density and ionization intensity due to the neutral gas diffusion motion after the ionization balance is broken.

Figure 7. Simulation results of thruster performance with step response of MW. (a) Screen grid current and (b) thrust.

Preliminary experiments have validated that thrust noise can be obviously reduced by feedback control when the MW power and ion beam current are respectively chosen as the controlled variable and feedback parameter (Zhou et al. Reference Zhou, Pang, Liu, Zhou and Niu2022). Therefore, four first-order inertia response processes of MW power, shown as (3.1)∼(3.4), are used to further study their influences on the SG current. Their response times are respectively 1 ms, 5 ms, 10 ms and 15 ms, which are chosen considering that the characteristic time of the SG current with the step variation of MW power is approximately 10 ms. To further study the influence of MW power overshoot characteristics on the SG current, the second-order under-damped processes with the same response time are used as the model input, shown as (3.5)∼(3.8)

(3.1)\begin{gather}{G_p} = \frac{{5000}}{{s + 5000}},\end{gather}

(3.2)\begin{gather}{G_p} = \frac{{1000}}{{s + 1000}},\end{gather}

(3.3)\begin{gather}{G_p} = \frac{{500}}{{s + 500}},\end{gather}

(3.4)\begin{gather}{G_p} = \frac{{333.33}}{{s + 333.33}},\end{gather}

(3.5)\begin{gather}{G_p} = \frac{{3.6 \times {{10}^7}}}{{{s^2} + 6000s + 3.6 \times {{10}^7}}},\end{gather}

(3.6)\begin{gather}{G_p} = \frac{{4 \times {{10}^6}}}{{{s^2} + 2000s + 4 \times {{10}^6}}},\end{gather}

(3.7)\begin{gather}{G_p} = \frac{{{{10}^6}}}{{{s^2} + 1000s + {{10}^6}}},\end{gather}

(3.8)\begin{gather}{G_p} = \frac{{4.444 \times {{10}^5}}}{{{s^2} + 666.7s + 4.444 \times {{10}^5}}}.\end{gather}

The SG current under various MW power response processes is shown in figure 8. It can be seen from figure 8(a) that the SG current shows overshoot characteristics although there is no overshoot for MW power under the first-order inertia process. The overshoot of the SG current declines as the MW power response process is prolonged. It can be concluded from the above analysis that the average atom density decreases due to the neutral gas diffusion after the ionization balance is broken by the rise of MW power. The characteristic time is approximately 10 ms. When the response time of MW power is less than 10 ms, the MW power changes little during the time interval from the response time of the MW to that of the thruster (10 ms). However, the neutral gas density and ionization intensity still fall in this time interval, which results in the decline of the SG current. It is the main reason why there is an overshoot for the SG current and its response time is almost unchanged under different MW response processes when the response time of MW power is less than 10 ms. Besides, the overshoot of the SG current declines as the MW power response process is prolonged because the continuous increase of MW power can weaken the decrease degree of the SG current caused by the fall of the average atom density. When the response time of MW power exceeds the characteristic time of the volume effect (10 ms), the increase of MW power can make the SG current increase continuously and therefore the overshoot phenomenon of the SG current gradually disappears.

Figure 8. Simulation results of screen current under different responses of MW. (a) First-order inertia process and (b) second-order under-damped process.

It can be seen from figure 8(b) that there are always overshoots for SG current under the second-order under-damped MW power response processes because the output MW power itself has overshoots. Under this condition, the response process of SG current shows similar characteristics to that under the first-order inertia process. When the response time of MW power is shorter than the characteristic time of the volume effect, there is a stage where the SG current still falls after the MW power becomes stable. This phenomenon gradually disappears when the response process of MW power is prolonged. Besides, it can be seen from figure 8(b) that the overshoots of the SG current decrease under slower MW power response processes, although the overshoots of MW powers are the same under these second-order under-damped response processes. This phenomenon reflects that the overshoots of the SG current depend on not only the overshoots of MW power, but also the response times of MW power, which are closely related to the characteristics of ECRITs.

3.2. Influence of flow rate on thruster dynamic characteristics

When the MW power and SG/AG voltage are set as 0.3 W and 500 V/−100 V, the SG currents under different flow rates is shown as figure 9, which shows that the simulation error is less than 5 %.

Figure 9. Comparison of experiments and model results with the variation of flow rate.

The state parameters in the discharge chamber are shown in figure 10 when the flow rate increases from 0.03 to 0.1 sccm in two steps. It can be seen from figure 10 that the average atom density increases with the step rise of the flow rate, which can prompt the ionization intensity. Therefore, electron density rises with the enhancement of ionization intensity, which consumes more electron energy and causes the drop of electron temperature. Besides, the electron density increment caused by the increase of flow rate from 0.03 to 0.06 sccm is more obvious than that from 0.06 to 0.1 sccm. Correspondingly, the atom density increment caused by the increase of flow rate from 0.03 to 0.06 sccm is less obvious than that from 0.06 to 0.1 sccm. It is mainly because the average electron temperature is higher under the first adjustment process, which causes that the enhancement of ionization intensity under the first adjustment process is more obvious than that under the second adjustment process. The response time of state parameters is 10 ms, consistent with the volume effect.

Figure 10. Simulation results of state parameters with step response of flow rate. (a) Atom density, (b) electron density, (c) electron temperature and (d) atom temperature.

It can be seen from figure 10(d) that the atom temperature firstly decreases and then rises after the step increase of the flow rate. The response time of the decrease process is 0.1 ms. It is mainly because the temperature of neutral gas flowing into the discharge chamber is below that previously in the discharge chamber. Accompanied by this process, the ion density increases with the enhancement of the ionization process. The collisions between ions and neutral gas heat the neutral gas and result in the increase of atom temperature. The response time of this process is approximately 10 ms. When the flow rate increases from 0.03 to 0.06 sccm, the ion density increment caused by the enhancement of ionization intensity is large, which leads to a strong heat effect of ions on atoms. Therefore, the atom temperature stays at a higher level after the flow rate increases. However, when the flow rate increases from 0.06 to 0.1 sccm, the heat effect of ions on atoms is weak, which causes the atom temperature to decrease due to inflow of low-temperature neutral gas. The response processes of the SG current and thrust with the step increase of flow rate are shown in figure 11.

Figure 11. Simulation results of thruster performance with step response of flow rate. (a) Screen grid current, (b) thrust and (c) sheath thickness.

It can be seen that the response characteristics of the thrust are similar to those of the SG current. It is mainly because the perveance condition of the thruster is under optimal perveance and the AG current is negligible compared with the SG current. Besides, the grid voltage is invariant during the adjustment process. Therefore, the thrust is mainly determined by the SG current. It can be seen that the response time of the SG current is approximately 10 ms, which is consistent with the characteristic time of the volume effect. Besides, it can be seen that the SG current decreases when the flow rate increases from 0.06 to 0.1 sccm. The SG current depends on the ion density upstream of the SG, ion Bohm velocity and ion extraction area. Although the increase of flow rate directly results in the rise of ion density upstream of the SG, the increment is small because the ionization intensity is limited by the low electron temperature. Meanwhile, the increase of flow rate results in the drop of electron temperature, which reduces the ion Bohm velocity upstream of the SG. Besides, it can be seen from figure 11(c) that the sheath thickness decreases with the increase of ion density upstream of the SG, which further results in the decline of the ion extraction area. Therefore, the SG current decreases when the flow rate rises from 0.06 to 0.1 sccm. When the flow rate increases from 0.03 to 0.06 sccm, the ion density increases at large scale. Therefore, the SG current still increases although the sheath thickness and Bohm velocity decrease under this adjustment process.

3.3. Influence of grid voltage on thruster dynamic characteristics

When the MW power and flow rate are set as 0.6 W and 0.06 sccm, the SG current under different SG voltages is shown as figure 12. The corresponding AG voltages are respectively −50 V, −50 V, −100 V, −150 V, −150 V and −200 V with the increase of SG voltages. It can be seen that the simulation results agree well with experiments, with an error of less than 5 %.

Figure 12. Comparison of experiments and model results with the variation of grid voltage.

The flow rate and MW power are respectively set as 0.06 sccm and 0.6 W. The state parameters in the discharge chamber are shown in figure 13 when the SG voltage steps from 200 to 500 V with 150 V intervals. The AG voltage remains unchanged during the first adjustment process and steps from −50 to −100 V during the second adjustment process. It can be seen that the electron density firstly decreases with the step increase of SG voltage because of the stronger extraction effect on the ions upstream of the SG by the grid system. The drop of electron density reduces the ionization intensity and less electron energy is consumed by the ionization process, which effects lead to the rise of the electron temperature. Ions in the discharge chamber move with the Bohm velocity, which is estimated to be 3 × 10³ m s⁻¹ if we assume the average electron temperature is 8 eV. The channel length of the thruster is 15 mm, which means that the characteristic time of this process is approximately 5 μs. The main drop of electron density occurs on a time scale that is of the order of this time interval. However, the decline of ionization intensity results in the increase of atom density. This response process is determined by the diffusion motion of neutral gas. Simulation results show that this process lasts approximately 0.5 ms. It can be seen that the response process of atom density due to step increase of the SG voltage is much faster than that due to the step increase of the MW power (figure 6c). It is mainly because the influence of grid voltages on state parameters in the discharge chamber is much weaker than that of MW power, which can be validated by the comparison between figures 6 and 13. Increase of atom density enhances the ionization intensity and results in the overshoot of electron density and electron temperature. Meanwhile, the response time of electron density and electron temperature is also prolonged to 0.5 ms because of the broken balance of ionization process and diffusion motion of neutral gas.

Figure 13. Simulation results of state parameters with step response of grid voltage. (a) Electron density, (b) electron temperature, (c) atom density, (d) atom temperature, (e) sheath thickness and (f) extraction area.

It can be seen from figure 13(d) that the variation characteristic of atom temperature is similar to that of ion density. It is mainly because the heat effect of ions on atoms by collisions dominates when the flow rate is stable. Besides, it can be seen the variation range of state parameters during the second adjustment process is larger than that during the first adjustment process, although the variations of SG voltages are the same. Analysis shows that it is mainly caused by the variation range of the extraction area. It can be seen from figure 13(e) that variations of sheath thickness are similar during the two adjustments. However, the sheath is thicker during the second adjustment process. Therefore, a similar variation of sheath thickness can result in a larger variation of the extraction area during the second adjustment, which can be validated from figure 13(f). Besides, it can be seen from figure 13(e) that the sheath thickness firstly increases and then falls, which forms an overshoot. The sheath thickness firstly increases mainly because of the stronger extraction capability and decline of ion density upstream of the SG with the step increase of grid voltages. Then the sheath thickness decreases mainly because the ion density upstream of the SG increases while the grid voltage remains invariant in the latter part of the response process.

The response processes of thruster performances are shown as figure 14. It can be seen that the response time of the SG current for the first and second adjustment processes are respectively 5 μs and 0.5 ms. The response time of the SG current during the second adjustment agrees well with that of state parameters, which is mainly influenced by the diffusion motion of neutral gas. However, the response time of the SG current during the first adjustment is consistent with the time interval when ions move through the discharge chamber with the Bohm velocity. The motion of neutral gas seems to have negligible influence on the response process of the SG current, which is caused by two reasons. On the one hand, the increments of state parameters are smaller under the first adjustment process compared with those under the second adjustment process. On the other hand, it can be seen from figure 13 that the variation direction of the ion density upstream of the SG is opposite to that of the ion Bohm velocity and ion extraction area. The SG current is determined by the product of these three parts. Therefore, variations of state parameters have negligible influence on the response of SG current when their ranges are inconspicuous.

Figure 14. Simulation results of thruster performance with step response of grid voltage. (a) Screen grid current, (b) AG current, (c) ion beam current and (d) thrust.

It can be seen from figure 14(b) that the AG current decreases to a very low level when the SG voltage increases from 200 to 500 V. When the SG and AG voltages are respectively 200 V and −50 V, the SG sheath is very thin and under the over-perveance condition. Under this condition, the ions ejected from the SG sheath tend to move into the sidewall of the AG and cause the AG current to be high. With an increase of grid voltages, the thruster gradually changes from over-perveance to optimal perveance and most ions ejected from the SG sheath can escape from the AG holes, which directly leads to the decrease of the AG current. The response time of the AG current for the first adjustment is approximately 5 μs, which is consistent with that of the SG current. There is an overshoot for the AG current, which is much higher than that of the SG current. The ion density in the grid system firstly decreases obviously with the increase of grid voltages. Then ions move from the discharge chamber to the grids, which contributes to the rise of ion density in the grid system and results in an overshoot of the AG current. There are no overshoots of the AG current for the second adjustment. This is mainly because most ions can escape from the AG holes when the SG sheath is under the optimal-perveance condition. The AG current mainly comes from the slow ions produced by CEX collisions downstream of the AG. Therefore, the ion density variation has negligible influence on the AG current for the second adjustment. It can be seen that the ion beam current and thrust have similar response processes with the step increase of grid voltages. The overshoot is obvious and its response time is approximately 5 μs for the first response process. The overshoot is inconspicuous and its response time is approximately 0.5 ms for the second response process. It is mainly determined by the response processes of the SG current and AG current.