## 3 Numerical Procedure  
The micronozzle operating regime depends on the ratio of the outlet pressure \(p_{2}\) to the inlet pressure \(p_{1}\). When \(p_{2}\) is equal to the ambient pressure and \(p_{2}/p_{1}\) is equal to or less than a critical pressure ratio \(\eta_{cr}\), the flow becomes sonic at the nozzle throat and it continues accelerate in the diverging section. The flow is fully expanded and a supersonic flow velocity forms at the nozzle exit. This is the normal operating regime of a supersonic nozzle and is desirable in many applications. Moreover, \(\eta_{cr}\) is a significant value for analyzing the flow in the nozzle as it changes from a subsonic to a supersonic velocity. \(\eta_{cr}\) can be calculated using the following expression from a one-dimensional isentropic relation:  
\[\eta_{cr}=\frac{p_{2}}{p_{1}}=\left(\frac{2}{k+1}\right)^{\frac{k}{k-1}}=0.487 \tag{14}\]  
The physical model and grid model in this paper are shown in Figure 1. It is a typical convergent-divergent nozzle and is based on the shape parameters of the nozzle used in the experimental study of [PERSON] [19] and the numerical study of [PERSON] et al. [3]. The angles of the converging and diverging parts are \(\alpha=45^{\circ}\) and \(\beta=20^{\circ}\), respectively. The nozzle throat half-width \(h\) is 100 \(\upmu\)m. Local grid refinement was carried out in the convergent part, divergent part and the nozzle throat. The simulation was carried out in the fully expanded condition and the flow was assumed to be steady and two-dimensional. The flow field of argon was driven by the pressure ratio of the outlet and inlet, and the inlet pressure was 101,325 Pa, with the initial outlet and ambient pressure of 3739 Pa. The wall temperature was set to 300 K. The average velocity at the throat was near the local speed of sound \(c^{*}=\sqrt{kR_{g}T}=337\) m/s, with the Reynolds number \(Re=\rho uh/\mu\) equaling 178. It was regarded as a laminar flow when the Reynolds number in the nozzle throat was less than 1000 [20]. The flow was assumed to be compressible and the density-based solver was employed. The convective terms were calculated based on a MUSCL (monotonic upstream-centered scheme for conservation laws) reconstruction, and the flow variable unit centered on the boundary of the grid was reconstructed with third-order accuracy [21]. The Roe-FDS was selected as the flux type. The double-precision floating point operation was chosen. Moreover, the convergence criterion for iteratively solving continuity, momentum and energy equations was set at \(10^{-8}\).  
The following assumptions were used for the discrete phase. On the one hand, it was assumed that the particles at the entrance were evenly distributed, and the type of nozzle wall was set to trap. This means that particles would be removed from the simulation if they hit the wall. On the other hand, one-way coupling was adopted owing to the volume fraction of particles being less than 10%; thus, the interaction between particles was ignored as well as the reverse influence of the particles on the fluid. The combined approach of computational fluid dynamics and the discrete particle method (CFD-DPM) was employed to investigate the motion characteristics of micro-/nano-particles in this study [22; 23]. When submicron particles collided with the nozzle wall, they were considered to stick to the surface. Particles are seeded with the gas velocity at the left boundary of the simulation domain. All numerical simulations were carried out with the commercial software package ANSYS FLUENT15.0.