Comparison of numerical predictions of the supersonic expansion inside micronozzles of micro-resistojets

. The present work provides a numerical investigation of the supersonic flow inside a planar micronozzle configuration under different gas rarefaction conditions. Two different propellants have been considered, namely water vapor and nitrogen, which relate to their use in VLMs (the former) and cold gas microthrusters (the latter), respectively. Furthermore, two different numerical approaches have been used due to the different gas rarefaction regime, i.e. the typical continuum Navier-Stokes with partial slip assumption at walls and the particle-based Direct Simulation Monte Carlo (DSMC) technique. As a result, under high-pressure operating conditions, both water and nitrogen flows supersonically expanded into the micronozzle without chocking in combination with a linear growth of the boundary layer on walls. However, when low-pressure operating condition are imposed and a molecular regime is established inside the micronozzle, a very rapid expansion occurred close to the nozzle exit in combination with a strong chocking of the flow and a micronozzle quality reduction of about 40%. Furthermore, water exhibited specific higher specific impulse than nitrogen above 60%.


Introduction
In the space field, the recent interest in the micro-power systems has driven the efforts of both scientific and industrial communities towards the development of small spacecrafts, thanks to the advancements in micro-fabrication technology of Micro Electro-Mechanical System (MEMS) devices. The design of such micro-and nano-satellites must consider the installation of a secondary propulsion system, which should produce small thrust forces, from few micronewtowns up to some millinewton, and high specific impulse, in compliance with severe mass, volume and power consumption constraints. In this concern, MEMS-based micro-propulsion systems, or microthrusters, represents a promising solution thanks to the possibility to be easily integrated on small satellites [1].
Among all micro-propulsion concepts, micro-resistojets represent an interesting choice because of the reduced system complexity and the simplicity of the working principle. Based on the propellant storage, they can be classified in vaporizing liquid microthruster (VLM) [2] and cold gas microthruster (CG). With respect to CGs, VLMs use liquid propellant, which can be easily stored and handled reducing weights and mass of the feeding system. However, some electrical power needs to be consumed for propellant vaporization and overheating. Based on the gas rarefaction condition established inside the micro-resistojet, these microthrusters can also be distinguished into i) high-pressure micro-resistojets (HPMR), and ii) low-pressure micro-resistojets (LPMR) [3], best known as free molecular micro-resistojets (FMMR) [4]. The choice between them is strictly related to the desired thrust level: the thrust level reduces up to few micronewtons when using FMMR, while HPMR can provide thrust forces of some millinewtons.
By considering the global performance of micro-resistojets, the micronozzle overall efficiency is strongly affected by the entity of the viscous effects [5,6]. These last ones are influenced by the degree of gas rarefaction, since it determines the mechanisms of interaction between gas-gas molecules and solid wall-gas molecules. Furthermore, the low-pressure condition into FMMR strengthens the gas rarefaction condition. In general, in micronozzles the continuum assumption and the no slip condition at walls are violated and two different regimes could establish based on the entity of the Knudsen number, which is defined as the ratio between the mean free molecular path and the characteristic length. In particular, when 0.01<Kn<0.1 the slip flow regime, while the transitional flow regime occurs if 0.1<Kn<10. In cases involving the first regime, the Navier-Stokes (NS) equations are still valid for numerical modeling, even though in combination with partial slip models at walls [7]. This condition is typical for HPMR. Instead, the transitional flow regime usually takes place into FMMR, and gas kinetic schemes such as Direct Simulation Monte Carlo (DSMC) [8] are required, as demonstrated in [9].
In this context, the present work aims to provide a numerical analysis of the performance of a planar tronco-conical micronozzle of MEMS-based micro-resistojets operating as VLM or CG, and both water vapor and gas nitrogen have been considered. the investigation involved both high-pressure and low-pressure operating conditions. Numerical 2D computations, using both NS with Maxwellian partial slip at walls with a tangential momentum accommodation coefficient (TMAC) equal to 0.8, and DSMC computations have been performed in 2D configuration depending on the gas rarefaction regime. global performance of the micronozzle in terms of thrust force, specific impulse, discharge coefficient and Isp-efficiency have been estimated. A qualitative description of the flow behaviour has been also provided in order to highlight differences and peculiarities in numerical solutions among the investigated cases.

Numerical approach
Numerical investigations were performed by using the open source CFD tool-box OpenFOAM© version 3.0.1. The geometry of the micronozzle was based on the VLM developed by Cen et al. [10]. In particular, the micronozzle is 120 μm depth, with an inlet width of 1070 μm, a throat width of 150 μm and an exit cross section of 1760 μm. The convergent angle is 45° while the divergent cone has 15° angle. Furthermore, a radius of curvature equal to 75 μm characterized the throat section, while a mixing region of 180 μm length was added before the entrance into the convergent region. At the inlet section, 8 microchannels with 90 μm width have been considered instead of the original 9 microchannels owning 80 μm width. Two different propellants have been considered, namely water vapor and gas nitrogen. In NS computations, the mass flow rate was imposed at the micronozzle entrance equal to 5.00x10 -6 kg/s for the H2O and 5.00x10 -6 kg/s for N2, in order to have an inlet pressure of about 2.2x10 5 Pa for both propellants, while in all DSMC computations the pressure of 300 Pa was directly set at the inlet with a mass flow rate equal to 4.14x10 -9 kg/s for the H2O and 7.27x10 -9 kg/s for N2.
The static temperature at the inlet and at walls was set to 300 K for gas nitrogen, while water vapor inlet and wall temperatures were set to 505.58 K, as retrieved by solving the steadystate boiling flow inside the heating chamber using an in-house 1D model [6]. The downstream pressure was set to 500 Pa in NS simulations, and 0.001 Pa for DSMC computations under low-pressure operating condition. Simulation were stopped after 2x10 -4 s since the steady state regime established into the micronozzle.

Navier-Stokes modeling
The density-based solver rhoCentralFoam [11] was used for compressible Navier-Stokes computations, in combination with the laminar turbulent model and Peng-Robinson equation of state [12] The central upwind scheme of Kurganov and Tadmor [13] was used for the flux terms and the Total Variation Diminishing (TVD) van Leer limiter [14] for interpolation. The time step was determined based on a maximum Courant number Comax of 0.2. The viscous governing equations were solved using a Preconditioned Conjugate Gradient/Diagonal Incomplete Cholesky scheme with a residual tolerance of 1 ×10 −8 . A partial slip boundary condition at walls was used since the maximum Knudsen number Kn into the divergent was strictly below 1×10 −1 . The tangential momentum accommodation coefficient (σTMAC) was set to 0.80 in combination with a first order Maxwell slip model, in reference of a generic gas flow on polished silicon microchannels. The computational domain is shown in Fig.1 (a). It extends 15Wexit downstream and 6.5Wexit upward, where Wexit is the width of the exit section.
The final mesh size is 23391 cells. The validation and verification of mesh and numerical modelling has been previously performed as described in [6] based on experimental data by Cen [10].

DSMC modeling
DSMC computations were carried out using the solver dsmcFoam+ [15]. The Variable Hard Sphere (VHS) model coupled with the Simplified Bernoulli Trials (SBT) model [16] defined the binary collisions between particles. Walls were treated as diffusive, and the Maxwellian thermal model determined the particle-wall collisions. The molecular properties used in DSMC computations are: the molecule diameter dmol (equal to 4.17x10 -10 for N2 and 5.78x10 -10 for H2O), the mass Mmol (equal to 4.65x10 -26 for N2 and 2.99x10 -26 for H2O), the temperature coefficient of viscosity ω (equal to 0.74 for N2 and 1 for H2O), the rotation DoF degrees of freedom (equal to 2 for N2 and 3 for H2O) and the rotation DoF degrees of freedom (equal to 1 for N2 and 3 for H2O). In DSMC computations, the number of simulated particles is usually less than the real number of particles. Consequently, an equivalent number of particles Np,eq need to be defined. In our simulations, we ensure the condition Np,eq > 10. Therefore, the total number of simulated particles is strictly related to the mesh size, and it plays an important role on the accuracy of the simulation and they are estimated according to [17]. In particular, the grid step must satisfy the condition ∆ ≤ is the free molecular path, kB is the Boltzmann constant, dmol is the molecule diameter, p and T are the operating pressure and temperature respectively. Instead, the time step must be less than the mean time between collisions as follows Δ = � �� � � �̅ , where ̅ = �8 � ��� ⁄ is the mean thermal velocity, and ξΔt is the time step parameter. Based on equations above, the mesh size results to be 47080 cells (see Fig.1 (b)), while the time step was set to 4x10 -9 s corresponding to the condition ξΔt ≤ 0.2, in accordance with typical DSMC requirements [18]. The methodology in DSMC computation has been previously validated based on experiments performed in [19] on a tronco-conical axisymmetric micronozzle.

Performance estimation
The thrust force F was estimated by means of a cell-based integration at the exit section, as a result of the sum of two terms, namely the jet and the pressure thrusts. The specific impulse was computed as �� = / � , where g0=9.81 m/s is the gravitational acceleration at sea level. The overall micronozzle efficiency � = � � , better known as micronozzle quality, is given by two different terms: the mass flow rate losses, measured by the discharge coefficient � , and the specific impulse losses represented by the I-sp efficiency � . In particular, the discharge coefficient was defined as � =̇̇� �� ⁄ , where the ideal mass flow rate ̇� �� was retrieved from the ideal rocket theory (IRT) [20] based on chocked nozzle condition and the isentropic flow hypothesis. Instead, the Isp-efficiency was estimated by imposing the momentum balance into the micronozzle [21] as follows: In cases involving water vapor, the thermal efficiency �� , defined as the ratio between the propulsive power and the thermal power required for vaporization and overheating, has been estimated as follows : where �,� and �,� are the heat capacities of the liquid and vapor phases, Ti,l is the feeding temperature of the liquid propellant, Tv is the saturation temperature and hL is the latent heat of vaporization. Figure 3 shows the contour of Mach number and temperature. Under high-pressure operating conditions, both water (test case 1) and nitrogen (test case 2) flows supersonically expanded into the micronozzle without chocking due to the boundary layer growth. This is mainly due to the 2D approximation, since the boundary layer growth along the lateral walls does not obstruct the flow passage, as confirmed by the linear growth of the displacement thicknesses retrieved from Eqs (2) and (3). By comparing water and nitrogen, results pointed out that water is more sensitive to the viscous boundary layer growth at the same inlet pressure condition. This is probably due to the higher inlet and wall temperatures, as confirmed by the stronger thermal boundary layer growth in Fig. 3 (b).

Results
When low-pressure operating condition are imposed and a molecular regime is established inside the micronozzle, the supersonic expansion of the flow is completely different, as shown in Fig. 5. In particular, the sonic flow condition is reached into the divergent region, towards the exit section, and the exit Mach number significantly reduces to a maximum of about 1.4. Consequently, the boundary layer growth is affected by the strong gas rarefaction condition: due to the very rapid expansion close to the nozzle exit, a strong chocking of the flow occurs as evinced by the momentum thickness estimated based on Eqs. (2) and (3). This negatively impacts on the mass flow rate losses, as it will be shown and discussed in the following section 3.4. Also, the thermal boundary layer growth is not more linear, and the heating provided by walls seems to play a crucial role in performance degradation as evinced by the different temperature contour plot into the divergent section in Figures 3 (b) and 4(b). In addition, it figures out that the difference in Mach and temperature fields between water (test case 3) and nitrogen (test case 4) diminishes. Table 3 summarizes the micronozzle performance in all investigated gas rarefaction conditions. By comparing propellants, water performs better than nitrogen under highpressure condition. In fact, the specific impulse of water (test case 1) is about 65% higher than nitrogen (test case 2). This is confirmed by the thrust force: in fact, a similar thrust force is produced but the propellant consumption is reduced of about 43% when using water. When a stronger gas rarefaction condition occurs into the micronozzle, the specific impulse reduced of about 31.5% and 29% for water (test case 3) and nitrogen (test case 4). Conversely, the Ispefficiency is not significantly affected when low-pressure operating condition are imposed.  4). Finally, it is worth to observe that ,when using water, the stronger the gas rarefied condition is, the lower is the power consumption due to the lowering of the mass flow rate coupled with severe losses in thermal efficiency.

Conclusions
The present work presents a numerical investigation of the performance of a planar micronozzle of MEMS-based micro-resistojets. Two fluids have been considered: water vapor in relation to a vaporizing liquid microthruster, and gas nitrogen, which refers to a cold gas microthruster. Furthermore, two different gas rarefaction conditions have been simulated, in reference to a high-pressure operating condition and a low-pressure operating condition. Due to the different gas rarefaction regimes, both NS and DSMC computations have been performed. The former to simulate the supersonic expansion at high-pressure condition/low gas rarefaction, the latter to simulate the supersonic expansion at low-pressure condition/transitional regime. Results highlighted that, under high-pressure operating conditions, both water and nitrogen flows supersonically expanded into the micronozzle without chocking in combination with a linear growth of the boundary layer on walls. Instead, when low-pressure operating condition are imposed and a molecular regime is established inside the micronozzle, the boundary layer growth is affected by the strong gas rarefaction condition: due to the very rapid expansion close to the nozzle exit, a strong chocking of the flow occurs as evinced by the momentum thickness. The overall performance analysis quantified the losses in terms of thrust, specific impulse and micronozzle quality. In particular, water performs better than nitrogen under high-pressure condition with a higher specific impulse of about 65%. When the stronger gas rarefaction condition occurs into the micronozzle, the specific impulse reduced of about 31.5% and 29% for water and nitrogen. Furthermore, the micronozzle quality significantly reduced mainly due to the mass flow rate losses represented by the discharge coefficient which dropped up to 40%.