Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 56, 3, pp. 713-725, Warsaw 2018 DOI: 10.15632/jtam-pl.56.3.713 STUDY ON THE EFFECTS OF PREHEATED WALL/PLATES IN MICROTHRUSTER SYSTEMS Elyas Lekzian, Hamid Parhizkar, Asghar Ebrahimi Aerospace Faculty, Malek Ashtar University of Technology, Tehran, Iran e-mail: lekzian@mut.ac.ir; hparhiz@mut.ac.ir; ebrahimi@mut.ac.ir In the present paper, effects of pre-heated walls/plates on microthrusters performance are studied using a DSMC/NS solver. Three microthruster configuration types are studied. Type 1 is a cold gas microthrster. Microthruster type 2 has pre-heated walls. Pre-heated plates are inserted inside the chamber of microthruster type 3. It is observed that in mi- crothruster type 2 the flow is accelerated and the specific impulse is elevated. However, by insertion of the pre-heated plates inmicrothruster type 3, viscous effects have stronger nega- tive influence and the thrust is decreased. By implementing temperature gradients on walls in type 2 and on plates in type 3, it is observed that a higher temperature gradient enhances performance parameters ofmicrothruters. Among all types ofmicrothrusters,microthruster type 2 with pre-heated walls has the highest thrust and specific impulse. Microthruster ty- pe 3with a temperature gradient of 300-500Khas theminimumthrust due to a considerable decrease in the mass flow rate. Keywords: DSMC/NS solver, pre-heated walls/plates, microthruster, performance parame- ters, temperature gradient 1. Introduction Lots ofmicro devices such asmicrochannels, micro heat sinks,microturbines,microengines, and microthrusters were developed with advancement in facbrication. Due to advantages of micro- -electro-mechanical-systems (MEMS) compared to theirmacro counterparts, they find increased applications in a variety of industrial andmedical fields (Gad-el-Hak, 2005a,b). Theseminiatu- rized devices utilize smaller volumes while offering the possibility of parallel operation on a chip andreduce the riskof thewhole systemfailure. Inaddition toadvantages of scaleminiaturization, they utilize little energy, offer high sensitivity andworkwith great accuracy. Besides the present applications of MEMS in electrical, structural, fluidic, transport and control aspects (Gad-el- -Hak, 2005a,b), the potential applications such as attitude control of small satellites usingmicro thrusters in deep space (Janson et al., 1999; Platt, 2002; Osiander et al., 2005; Mihailovic et al., 2011) are a field of interest. In 1960s, low thrust devices were analyzed experimentally (Milligan, 1964; Rothe, 1971). Then numerical methods were proposed by Rae (1971) and Boyd et al. (1992) in order to inve- stigate converging-diverging nozzles numerically. Boyd used a DSMC method and proved that this method was accurate enough to analyze the flow especially in the expanding section of the thruster with a near vacuum exit boundary conditions. Small nozzles were studied by Bayt (1999). Hemodeled aMEMSbased nozzle by usingNavier-Stokes simulations and observed that by a decrease of Reynolds number the propulsive efficiency was decreased due to an increase in boundary layer thickness.Amicrothruster utilized inGP-B spacecraftwas studiedbyJafryt and Beukelt (1994). They studied back pressure effects on the behavior of flow in the nozzle. Ivanov et al. (1999) studied the effects of the throat Reynolds number on the specific impulse. They demonstrated that an increase of losses could led to over-prediction of the specific impulse at 714 E. Lekzian et al. the exit of the nozzle. Alexeenko et al. (2002) investigated an axisymmetric 3D micronozzle by theDSMCmethod andNavier-Stokes solver. Their results showed that gas expansion increased the specific impulse. Xie (2007) demonstrated that flow in MEMS nozzles was simulated accu- rately by the DSMCmethod and the Navier-Stokes equations even when adding slip boundary conditions became invalid when the average Knudsen numberwas about 0.01 in his special case. Surface roughness of themicro nozzle was studied by Torre et al. (2010). They observed shocks near thewalls due to roughness. Sun et al. (2009) conducted aDSMC-FVMmethod to simulate the flow inside amicro nozzle. They investigated the effect of inlet pressure on the flowfield insi- de the nozzle. They concluded that when the inlet pressurewas increased, the distance between the throat and the point that propellant velocity surpassed sonic velocity became smaller. In the present study, pressure driven microthrusters are investigated. Numerical simulation of the flow passing through the thruster is introduced first. ADSMC/NS solver is utilized. Due to lack of experimental data atmicroscales, the solver is first verifiedby a simulatedmicronozzle. The solver is then utilized to simulate microthrusters. Effects of heating the flow by increasing thewall temperature and by insertion of high temperature plates inside the domain are studied, and performance parameters are then investigated. Next, the effects of the temperature gradient on the performance parameters are analyzed. At the end, the thrust and specific impulse of all microthrusters configurations are compared. 2. Basic theory 2.1. DSMC solver DSMCmethod enables numerical solution of theBoltzmann equations. Themethod benefits from the statistical solution of the particles behavior represented by kinetic theory (Bird, 1994). For an accurate DSMC simulation, four issues must be considered: The cell size, time step, number of particles per cell and the mean collision separation distance. The cell size must be 1/3 of the mean free path (Hadjiconstantinou, 2000; Pfeiffer et al., 2013) to prevent errors in diffusion. In the present study, the method proposed by Nance et al. (1998) is utilized to determine the grid distribution. The time step must be smaller than the mean collision time in order to properly distinguish the particle freemovement from a pair collision. In order to obtain physical and meaningful solutions, Bird (2007) propose that the particle per cell (PPC) would be around 7.However, Le et al. (2006) showed that for the low speed portion of the flownear the inlet of the computational domain, an average of 20 particles per cell gives precise solutions. In this paper, the study is begunwith PPC=7, and then a PPC independence study is conducted. Another parameter that affects the results is the mean collision separation distance. According to Moss and Bird (2005), the mean separation distance has to be smaller than the mean free path in each cell throughout the computational domain. The DSMFOAM solver of OPENFOAM-2.3.1 is used in DSMC simulations. The boundary condition proposed by Liu et al. (2006) is applied in simulation. The boundary condition in microchannel systems strongly affects the results, as described by Nance et al. (1998). Current boundary treatment is simple and stable. Many publications (Liou and Fang, 2000; Wang and Li, 2004; Roohi et al., 2009) demonstrated the accuracy of this boundary condition. The solver is capable of parallel operation and can model any geometry and any number of gas species. For calculation of post collision velocities andmodeling of particle collisions, the Variable Hard Sphere (VHS) model is used. This model is widely applied because of its simplicity and good approximation of intermolecular collisions (Le et al., 2006). Study on the effects of preheated wall/plates in microthruster systems 715 2.2. NS solver By using the Chapman-Enskog procedure, the Navier stokes equations can be derived from the Boltzmann equations. The conservation of mass, momentum and energy are used (OPENFOAM, 2014).Maxwell velocity slip boundary conditions including thermal creep effects are implemented on the walls and plates (O’Hare et al., 2007) Vg−Vw =− 2−σu σu λ µ τ − 3 4 Pr(γ−1) γp q (2.1) where σu is the momentum accommodation coefficient. Subscripts g and w denote the gas adjacent to the wall and the wall, respectively. The parameter γ is the specific heat ratio and Pr is the Prandtl number. The shear stress is τ =S ·(n ·Π) and the heat flux is q=Q ·S. The tensor S= I−nn, where I is the identity tensor, and n is the normal unit vector to the wall. Q is the heat flux vector along the wall, andΠ is the stress tensor at thewall. The parameter µ is the gas viscosity and can be calculated using both the power law and the Sutherland law. The power law model is slightly less accurate than the Sutherland law at moderate temperatures. Therefore, the Sutherland law is used for calculation of gas viscosity as follows (Le and Roohi, 2015) µ=AS √ T3 T +TS (2.2) The parameter AS and TS are constants, AS = 1.41 · 10 −6Pas/K2, TS = 111K for nitrogen. Maxwellian definition can be used for calculation of the molecular mean free path λ (Le et al., 2012), but for theVHS collisionmodel themolecularmean free path can be calculated as follows (Bird, 1994) λ= 2(5−2ω)(7−2ω) 15 √ m 2πkT µ ρ (2.3) whereω is the temperature exponent of the viscosity coefficient in thepower lawviscositymodel, ω = 0.74 for nitrogen (Bird, 1994). The parameter m is molecular mass, and k is Boltzmann constant. It has been observed that the temperature of the rarefied gas adjacent to the wall is not equal to the wall temperature. Therefore, temperature jump boundary conditions are implemented on the walls when a dilute gas is simulated by the NS equations. In the present study, the second order Smoluchowski temperature jump boundary condition is used as follows (Karniadakis et al., 2006) Tg −Tw =− 2−σT σT 2γ γ+1 1 Pr λ ∂T ∂n − 2−σT σT 2γ γ+1 1 Pr λ2 2 ∂2T ∂n2 (2.4) where σT is the thermal accommodation coefficient and n is the unit vector normal to the wall. Momentum and thermal accommodation coefficients are equal to unity in order to simulate diffuse reflector walls. For NS simulation of the flow, RhoCentralFoam solver in OPENFOAM-2.3.1 is used. RhoCentralFoam is a density-based compressible flow solver based on central-upwind schemes of Kurganov and Tadmor. 3. Solver algorithm Tomakeuseof theDSMC/NSsolver, aparameter is required tobedefined in order todistinguish the continuumandmolecular regions. In the present study, the localKnudsen number is utilized to distinguish rarefied and continuum regions. The local Knudsen number is defined as below 716 E. Lekzian et al. Knlocal = λlocal ϕ ( ∆ϕ ∆x ) −1 (3.1) where ϕ is one of the macroscopic properties in each cell including: velocity, temperature, or density. In the NS solver, the parameter λlocal is calculated using equation (2.3). In the DSMC solver, the local mean free path is determined as below (Bird, 1994) λlocal = 1 √ 2πd2n (3.2) where d = 4.17 · 10−10m is the nitrogen molecular diameter, and n is the density num- ber within each cell. The local Knudsen number is determined as follows: Knlocal = max(KnlocalT ,Knlocalu,Knlocalρ) (Boyd et al., 1995). Based on the procedure proposed by Schwartzentruber et al. (2007), the DSMC/NS solver simulation cycle is as follows: • Entire computational domain is solved by the NS solver, • Local Knudsen number is calculated and the continuum and rarefied regions are distingu- ished, • Rarefied regions are solved by the DSMC solver, • Local Knudsen number is recalculated, and the interface of the continuum and rarefied region is changed, • Steps 3 and 4 are repeated until location of the interface region does not change. 4. Validation 4.1. DSMC solver validation There is no experimental study onmicrothruster systems atmicro scales. Therefore, in order to verify the accuracy of the DSMC solver, simulation results are compared with Liu et al. (2006). After a careful grid study based on the method proposed by Nance et al. (1998), the domain is divided into 400 cells in the x direction and 150 cells in the y direction by using a structured grid. The temperature and Mach number contours are depicted in Fig. 1a. The centerline temperature andMach number are presented in Fig. 1b, and the results are compared with those presented in Liu et al. (2006). Considering that the mean separation distance between the simulation particles must be smaller than the mean free path in each cell, two different particles per cell (PPC) (PPC=15, PPC=20) are chosen. As seen in Fig. 1b, a similarMach number and temperature distribution is achieved. Therefore, correct simulation results are obtained bybothPPC=15andPPC=20. Also, agreements of the results with those of Liu et al. (2006) shows that the DSMC solver is accurate and the results are trustworthy. 4.2. DSMC/NS solver validation In the previous Section, the entire domain is solved by the DSMC solver. In this Section, the DSMC/NS solver is used to obtain the solution. As mentioned earlier, the local Knudsen number needs to be determined. The domain is divided into 400×150 cells using a structured grid (Fig. 2a). Figure 2b shows the density based local Knudsen number which is themaximum local Knudsen number. It is observed that the boundary layer region in the divergent section of the micronozzle is in the rarefied regime since the local Knudsen number in this part is more than 0.05 (Boyd et al., 1995). Due to expansion of flow in the divergent section, the rarefied region is increased toward the end of nozzle. It is demonstrated by Darbandi and Roohi (2011) that solving the flow using the NS solver in the divergent section is not accurate even in the Study on the effects of preheated wall/plates in microthruster systems 717 Fig. 1. (a) Temperature, Mach contours and graphs of the micronozzle; (b) current temperature and Mach number compared with Liu et al. (2006) centerline. Therefore, the DSMC solver is used to simulate the entire divergent section and the NS solver is used for simulation of the convergent section of the micronozzle. The convergent section of the domain is divided into 100×150 cells and the divergent section is divided into 300×150 cells. Careful grid study is carried out based on the reference (Nance et al., 1998) to choose pre-mentioned grid division of the domain. Pressure and temperature are considered as inlet boundary conditions of themicronozzle at the convergent section. The solution at the end of the convergent section is considered as the inlet boundary conditions of the divergent section. This data is considered as the initial data for calculation of pressure, temperature and flux by averaging over the inlet of the divergent section. Fig. 2. Grid andKnudsen number distribution of the micronozzle: (a) structured grid, (b) Knudsen number contours using the DSMC solver (up) and the NS solver (bottom) Figure 3 shows temperature distribution along the micronozzle centerline achieved by the DSMC and DSMC/NS solvers. 718 E. Lekzian et al. Fig. 3. Temperature distribution using the DSMC and DSMC/NS It is observed that the result of the DSMC/NS solver is exactly the same as of the DSMC solver. Although the domain can be solved by the DSMC solver, the DSMC/NS solver is deve- loped to reduce computational cost. DSMC simulation of the abovemicronozzle (using an Intel Core i5 computer with 3GB RAM) takes 35 hours while using the DSMC/NS solver the time of computations is reduced to 23 hours (21 hours for the DSMC solver and 2 hours for the NS solver). 5. Microthruster problem statement In this Section, the DSMC/NS solver is used to simulate microthruster systems. The chosen microthruster device is shown in Fig. 4. The thruster consists of a rectangular channel which is connected to a converging-diverging nozzle. A buffer zone is attached to the end of the diverging nozzle to consider variation of parameters at the nozzle exit (Wu et al., 2001). Three types ofmicrothrusters are considered. In type 1, a cold gasmicrothruster is simulated (Fig. 4 without heater plates). Type 2 is the case when the flow is heated by increasing the wall temperature up to 900K (see Fig. 4 without heater plates, the wall temperature is 900K). The increasing of the wall temperature can be performed by using hot wire coated walls (Kundu et al., 2013). For type 3, heater plates (with temperature 900K) are inserted into the thruster (see Fig. 4 with heater plates). The overall geometry (but not the exact dimensions), the structure and idea of heating the fluidflow in type 3 is adapted fromKundu et al. (2013), Hitt et al. (2001) but not exactly with the same components. They fabricated and tested such a microthruster. Geometrical details of the three types of microthrusters are mentioned in Table 1. Table 1.Dimensional parameters of the three types of microthrusters L1 L2/L1 L3/L1 L4/L1 L5/L1 H1 H2/H3 H3/H1 Type 1 – Type 2 0.6mm 0.4 0.62 0.83 0.17 0.3mm 0.5 0.33 Type 3 0.17 The gas inside the thrusters is nitrogen. In all types of the thrusters mentioned above, the inlet temperature of 300K is imposed. Since the velocity at the thruster exit is supersonic, no outlet pressure is set at the exit. Study on the effects of preheated wall/plates in microthruster systems 719 Fig. 4. Geometry of the microthruster 5.1. Flow properties comparison of type 1, 2, and 3 Temperature contours for the three types of microthrusters are shown in Fig. 5. Type 1 is considered as the reference case. In type 1, the temperature is decreased in the divergent section due to gas expansion. In type 2, the temperature is increased and reaches to the maximum value (less than 900K) in the centerline and then it decreases. In type 3, the temperature is increased (to the maximum temperature of 900K of pre-heated plates), it remains constant in the pre-heated plates section and then is decreased. Fig. 5. Temperature contours for three types of microthrusters: (a) type 1 – cold gas, (b) type 2 – pre-heated walls, (c) type 3 – pre-heated plates For accurate analysis of types 1, 2, and 3, the flow properties along the microthruster cen- terline are presented and compared in Fig. 6. These graphs illustrate the effect of pre-heated walls/plates inside the device on the gas flow properties. A significant pressure reduction occurs when the flow passes over heater plates in type 3 while the pressure along the centerline chan- ges more smoothly in type 2 and 1. At the beginning and at the end of the heater plates the temperature is increased and decreased rapidly in type 3, whereas the temperature is increased more smoothly in type 2. In the convergent section, it is seen that the heater plates of type 3 microthruster systemaffect the downstream temperatureflowfield and increase the temperature downstream the plates. It can be concluded that the pressure loss at the heater plates in type 3 is converted into a temperature increase. 720 E. Lekzian et al. Fig. 6. Flow properties comparison of three types of microthrusters The temperature at the diverging section of the micro propulsion device type 2 is almost similar to the cold gas microthruster. The exit velocity for type 2 and 3 is identical and greater than type 1. Therefore, the heating process accelerates the flow. Heating the flow also affects theMach number at the exit. In type 2 and 3, the choking place of flow is not changed, but the exit Mach number of type 2 is greater than type 3. Such behavior is expected because the flow velocity at the exit is almost identical for type 2 and 3, but the temperature at the nozzle exit is higher in type 3 than type 2. Therefore, the exit Mach number of type 2 is greater than type 3. 5.2. Microthrusters performance results Thrust, specific impulse, exit velocity and the mass flow rate are compared for three types of microthrusters. The total thrust is the sum of thrusts of each cell at the exit. The thrust is calculated using F = ṁVexit. The specific impulse is calculated as Isp =F/(ṁg). Quick comparison of the performance of three preceding microthruster systems is provided in Table 2. It is observed that the thrust and the specific impulse are highest in type 2. The- refore, a higher temperature would result in a higher thrust and specific impulse. Meanwhile, the minimum of thrust occurs in type 3, where the heater plates are inserted in. In type 3, the viscous effects of heater plates dominate the gas expansion and reduce the thrust due to signifi- cant viscous effects of the plates. In type 2, the flow is heated and accelerated but there are no Study on the effects of preheated wall/plates in microthruster systems 721 heater plates in the domain to increase the frictional forces. Hence, the performance parameters of type 2 are higher than the two others. Table 2.Comparison of performance of three types of microthrusters ṁ [g/s] Vexit [m/s] F [mN] Isp [s] Type 1 12.3 980 11.1 106.7 Type 2 11.8 1310 14.7 144.7 Type 3 7.4 1290 7.1 121.3 5.3. Effects of temperature gradient Kundu et al. (2013) showed that there was a temperature gradient in a hydrogen proxide monopropellant microthruster at the chamber section. Therefore, in this Section, the effect of temperature gradient on the thruster performance is studied. Three cases of temperature gradients are implemented on thewall ofmicropropulsion device type 2 and on the heater plates of microthruster type 3 (see Table 3). Table 3.Comparison of performance of three types of microthrusters Type of Gradient case Temperature at the beginning Temperature at end of microthruster of wall/heater plates [K] wall/heater plates [K] Linear: 1 300 500 Type 2/type 3 Linear: 2 300 700 Linear: 3 300 900 Grid/PPC studies have been carefully carried out for all gradient cases. In this Section, only the results obtained from simulation are presented. Temperature gradients of microthrusters 2 and 3 without the gradient and with the gradient (case 3) are depicted in Fig. 7. Fig. 7. Comparison of temperature contours for gradient case 3 and constant temperature: (a) gradient case (up), constant pre-heated wall (900K) (bottom) for microthruster 2, (b) gradient case (up), constant pre-heated plate (900K) (bottom) for microthruster 2 It is observed that the exit temperature of microthruster 2 is slightly larger than the exit temperature ofmicrothruster 2 with gradient case 3. The same temperature behavior is seen for thruster type 3. Table 4 provides a quick review of the microthrusters performance. It is demonstrated that temperature gradient case 3 has a higher thrust and specific impulse for both types of microthrusters (type 2 and 3). 722 E. Lekzian et al. Table 4.Performance comparison of three types of microthrusters Type of Gradient case Mass flow Exit velocity Thrust Specific microthruster rate [gs−1] [ms−1] [mN] impulse [s] Linear: 1 12.1 1030 12.3 119.7 2 Linear: 2 12 1211 13.7 130.4 Linear: 3 11.7 1305 14.6 142.3 Linear: 1 6.85 1003 4.41 85.1 3 Linear: 2 6.9 1111 5.4 99.8 Linear: 3 7.2 1263 6.8 118.5 6. Performance comparison with/without temperature gradient Figure 8 compares the thrust, specific impulse, mass flow rate and the exit velocity of all types of microthrsueters with the constant temperature and with temperature gradient cases. Fig. 8. Comparison of performance parameters of three types of microthrusters with constant pre-heated walls/plates and with gradient cases It is observed that the highest thrust and specific impulse occurs inmicrothruster 2with the constantwall temperature of 900K. Insertion of pre-heated plates intomicrothruster 3 decreases the thrust.Pre-heated plates/walls increase the exit velocity. However, by implementing ahigher temperature gradient to walls/plates, the microthrusters performance is elevated. Study on the effects of preheated wall/plates in microthruster systems 723 7. Conclusion ADSMC/NS solver has been used in this paper. Accuracy of simulations has been first verified by simulation of a micronozzle. Then the flow inside three different microthruster system con- figurations has been simulated. A cold gas microthruster as type 1, pre-heated wall propulsion device as type 2, and a thruster with pre-heat plates inside the domain as type 3 have been simulated. It is observed that an increase in the fluid temperature accelerates the flow. Insertion of pre-heated plates inside the domain increases pressure loss, whereas heating up the walls of the thruster does not increase the frictional forces and, consequently, does not increase the pressure loss. It is demonstrated that type 2 has higher performance parameters than type 1 and 3. Also, temperature gradients have been implemented into the pre-heatedwalls and plates. It is observed that a higher gradient would elevate the thrust and specific impulse. Among all simulated cases, microthruster of type 2 with constant pre-heated walls (no temperature gra- dient) has better performance than all other cases. Microthruster 2 benefits from heating up which accelerates the flow, whereas it does not suffer from the frictional forces caused by the heater plates in type 3. References 1. Alexeenko A., LevinD., Gimelshein S., Collins R.,MarkelovG., 2002,Numerical simu- lation of high-temperature gas flows in a millimeter-scale thruster, Journal of Thermophysics and Heat Transfer, 16, 1, 10-16 2. Bayt R.L., 1999, Analysis, fabrication and testing of a MEMS-based micropropulsion system, PhD, Aerospace Computational Design Laboratory,Department of Aeronautics andAstronautics, Massachusetts Institute of Technology 3. Bird G., 2007, Sophisticated DSMC, Notes Prepared for a Short Course at the DSMC07 Meeting, Santa Fe, USA 4. Bird G.A., 1994,Molecular Gas Dynamics and the Direct Simulation of Gas Flows, USA, Oxford Science Publications 5. Boyd I.D., Chen G., Candler G.V., 1995, Predicting failure of the continuum fluid equations in transitional hypersonic flows,Physics of Fluids, 7, 1, 210-219 6. Boyd I.D., Penko P.F., Meissner D.L., DeWitt K.J., 1992, Experimental and numerical investigations of low-density nozzle and plume flows of nitrogen,AIAA Journal, 30, 10, 2453-2461 7. Darbandi M., Roohi E., 2011, Study of subsonic-supersonic gas flow through micro/nanoscale nozzles using unstructured DSMC solver,Microfluidics and Nanofluidics, 10, 2, 321-335 8. Gad-el-Hak M., 2005a,MEMS: Applications, USA, CRCPress 9. Gad-el-Hak M., 2005b,MEMS: Introduction and Fundamentals, USA, CRCPress 10. HadjiconstantinouN.G., 2000,Analysis of discretization in the direct simulationMonteCarlo, Physics of Fluids, 12, 10, 2634-2638 11. Hitt D.L., Zakrzwski C.M., Thomas M.A., 2001,MEMS-based satellite micropropulsion via catalyzed hydrogen peroxide decomposition, Smart Materials and Structures, 10, 6, 1163 12. Ivanov M., Markelov G., Ketsdever A., Wadsworth D., 1999, Numerical study of cold gasmicronozzle flows, 37th Aerospace Sciences Meeting and Exhibit, Reno, NV, U.S.A., American Institute of Aeronautics and Astronautics 13. Jafryt Y., Beukelt J.V., 1994, Investigation of nozzle and plume expansions of a small helium thruster,Rarefied Gas Dynamics: Space Science and Engineering, 160, 136 724 E. Lekzian et al. 14. Janson S.W., HelvajianH., Hansen W.W., Lodmell J., 1999,Microthrusters for nanosatel- lites, Second International Conference on IntegratedMicro Nanotechnology for Space Applications, Pasadena 15. KarniadakisG.E., BeskokA.,AluruN., 2006,Microflows andNanoflows: Fundamentals and Simulation, Springer NewYork 16. Kundu P., Sinha A.K., Bhattacharyya T.K., Das S., 2013, Nanowire embedded hydro- gen peroxide monopropellant MEMS thruster, Journal of Microelectromechanical Systems, 22, 2, 406-417 17. Le M., Hassan I., Esmail N., 2006, DSMC simulation of subsonic flows in parallel and series microchannels, Journal of Fluids Engineering, 128, 6, 1153-1163 18. Le N.T., Roohi E., 2015,A new form of the second-order temperature jump boundary condition for the low-speed nanoscale and hypersonic rarefied gas flow simulations, International Journal of Thermal Sciences, 98, 51-59 19. Le N.T., White C., Reese J.M., Myong R.S., 2012, Langmuir-Maxwell and Langmuir- -Smoluchowski boundary conditions for thermal gas flow simulations in hypersonic aerodynamics, International Journal of Heat and Mass Transfer, 55, 19, 5032-5043 20. LiouW., FangY., 2000, Implicit boundary conditions for direct simulationMonteCarlomethod inMEMSflowpredictions,CMES –ComputerModeling in Engineering and Sciences,1, 4, 119-128 21. Liu M., Zhang X., Zhang G., Chen Y., 2006, Study on micronozzle flow and propulsion per- formance using DSMC and continuummethods,Acta Mechanica Sinica, 22, 5, 409-416 22. Mihailovic M., Mathew T., Creemer J., Zandbergen B., Sarro P., 2011,MEMS silicon- -based resistojet micro-thruster for attitude control of nano-satellites, Solid-State Sensors, Actu- ators and Microsystems Conference (TRANSDUCERS), Beijing, IEEE: 262-265 23. MilliganM.W., 1964,Nozzle characteristics in the transition regimebetween continuumand free molecular flow,AIAA Journal, 2, 6, 1088-1092 24. Moss J.N., BirdG.A., 2005,Direct simulationMonteCarlo simulations of hypersonic flowswith shock interactions,AIAA Journal, 43, 12, 2565-2573 25. Nance R.P., Hash D.B., Hassan H., 1998, Role of boundary conditions in Monte Carlo si- mulation of microelectromechanical systems, Journal of Thermophysics and Heat Transfer, 12, 3, 447-449 26. O’Hare L., Lockerby D.A., Reese J.M., Emerson D.R., 2007, Near-wall effects in rarefied gasmicro-flows: somemodern hydrodynamic approaches, International Journal of Heat and Fluid Flow, 28, 1, 37-43 27. OPENFOAM, 2014),OPENFOAM: The Open Source CFD Toolbox 28. OsianderR., DarrinM.A.G.,Champion J.L., 2005,MEMS andMicrostructures in Aerospace Applications, USA, CRCPress 29. Pfeiffer M., Mirza A., Fasoulas S., 2013, A grid-independent particle pairing strategy for DSMC, Journal of Computational Physics, 246, 28-36 30. Platt D., 2002, A monopropellant milli-Newton thruster system for attitude control of nanosa- tellites,AIAA/USU Conference on Small Stellites, Utah 31. RaeW.J., 1971, Some numerical results on viscous low-density nozzle flows in the slender-channel approximation,AIAA Journal, 9, 5, 811-820 32. RoohiE.,DarbandiM.,MirjaliliV., 2009,Direct simulationMonteCarlo solutionof subsonic flow throughmicro/nanoscale channels, Journal of Heat Transfer, 131, 9, 092402 33. Rothe D.E., 1971, Electron-beam studies of viscous flow in supersonic nozzles,AIAA Journal, 9, 5, 804-811 Study on the effects of preheated wall/plates in microthruster systems 725 34. Schwartzentruber T.E., Scalabrin L.C., Boyd I.D., 2007, A modular particlecontinuum numerical method for hypersonic non-equilibrium gas flows, Journal of Computational Physics, 225, 1, 1159-1174 35. SunZ.-X., Li Z.-Y.,HeY.-L.,TaoW.-Q., 2009,Coupled solid (FVM)-fluid (DSMC) simulation of micro-nozzle with unstructured-grid,Microfluidics and Nanofluidics, 7, 5, 621-631 36. SunZ.-X.,TangZ., HeY.-L., TaoW.-Q., 2011,Proper cell dimension and number of particles per cell for DSMC, Computers and Fluids, 50, 1, 1-9 37. Torre F.L., Kenjeres S., KleijnC.R.,Moerel J.-L.P., 2010,Effects of wavy surface rough- ness on the performance of micronozzles, Journal of Propulsion and Power, 26, 4, 655-662 38. Wang M., Li Z., 2004, Simulations for gas flows in microgeometries using the direct simulation Monte Carlo method, International Journal of Heat and Fluid Flow, 25, 6, 975-985 39. Wu J., Lee W., Lee F., Wong S., 2001, Pressure boundary treatment in internal gas flows at subsonic speed using the DSMCmethod,Rarefied Gas Dynamics: 22nd International Symposium, USA, American Institute of Physics, 408-416 40. Xie C., 2007, Characteristics of micronozzle gas flows,Physics of Fluids, 19, 3, 037102 Manuscript received December 31, 2016; accepted for print November 7, 2017