Numerical simulation of parallel-plate particle separator for estimation of charge distribution of PM2.5

Abstract A numerical simulation of an instrument that is used to measure the charging state of PM2.5 is conducted in order to clarify its measurement uncertainty and to improve its performance. The instrument, a parallel-plate particle separator (PPPS), is designed to classify aerosol particles according to their charging states and measure their quantities. The trajectories of submicron particles in the PPPS are numerically analyzed using the Lagrangian particle tracking method, taking into account the Brownian force and the electrostatic force. First, it is confirmed that the deterioration in the classification accuracy observed in the experiment is due to Brownian diffusion. The optimal condition that improves the accuracy is investigated through a parametric study by varying the balance of flow rates at the inlets, the geometry of the inlet and exit sections, and the applied voltage. It is found that decreasing the flow rate of the central inlet for aerosol or narrowing the central inlet improves the accuracy. The dependence of the accuracy on the flow rate is found to be in accordance with the experimental results. For charged particles, an optimum voltage that maximizes the classification accuracy is found. On the basis of the simulation results, we propose a method to determine the charge distribution of aerosol from the number of particles counted at each exit of the PPPS. In the test assuming aerosol in the air, the charge distribution determined from the number count at the exits is found to perfectly agree with the charge distribution specified at the inlet. Copyright © 2019 American Association for Aerosol Research


Introduction
Particulate matter with an aerodynamic diameter of 2:5 lm or less in the air, called PM2.5, is considered to have adverse effects on human health by deposition on their respiratory organs (Takahashi 2003). To understand and suppress these adverse effects, many researchers have studied particle deposition on human organs, as reviewed, for example, by Kim et al. (2015). In respect of the charging states of particles, Cohen et al. (1998) found that singly charged particles with a diameter of 125 nm deposited on a human lung model five times more than non-charged particles. Ali, Reddy, and Mazumder (2008) also showed that charged particles in the diameter range 0.5-10 lm deposited on a replica of the human throat more frequently than non-charged particles. These tendencies have also been confirmed by numerical simulations. Majid et al. (2012) computed the deposition of charged particles in a stochastic lung model and confirmed that the deposition rate can be more than doubled by the inhalation of charged particles. In addition, they found that deposition is enhanced as the number of charges increases. Koullapis et al. (2016) showed, by numerical simulation with a realistic computational domain of a human throat from lungs, that electrostatic charge increases particle deposition up to seven times higher at the division from mouth to throat. Hence, information on the charge distribution of particles in a realistic aerosol is required in order to understand particle deposition on human organs in more detail.
It is well known that the charge distribution of aerosol particles under a stationary condition, that is, long-duration exposure of particles to numerous bipolar ions, follows Fuchs' stationary charge distribution (Fuchs 1963). In particular, Fuchs (1963) showed that the charge distribution of particles with a diameter greater than 0:6 lm matched the Boltzmann charge distribution. The theory of Fuchs (1963) has been confirmed by experimental observations. Liu and Pui (1974) showed that the measured charge distribution of 0.55-1:17 lm particles agreed with the Boltzmann charge distribution. Wiedensohler and Fissan (1991) and Alonso et al. (1997) confirmed that 2.5-100 nm particles followed Fuchs' stationary charge distribution. However, by using an integral mobility analyzer based on the theory of Knutson and Whitby (1975), Forsyth, Liu, and Romay (1998) showed that the charge distributions of particles were different from the Boltzmann charge distribution immediately after their generation in the laboratory. In addition, for nano-sized particles, Ahn et al. (2001) observed that the ratio of singly charged to doubly charged SiO 2 particles is higher than that from Fuchs' stationary charge distribution at particle diameters up to 50 nm in an H 2 /O 2 /TEOS flame. They used a method that employs a pair of differential mobility analyzers (DMAs) which can yield particles with a narrow range of electrical mobility (Hinds 1982), called the tandem DMA technique (Rader and McMurry 1986). This technique was refined by Kim et al. (2005) and Wild, Meyer, and Kasper (2012) for more accurate measurement. Further, some other methods and devices have also been developed to measure the charge distribution of particles (Kousaka et al. 1981;Buckley, Wright, and Henshaw 2008;Marra et al. 2009).
Recently, attempts have been made to measure the net charging states of ambient aerosols. Jayaratne, Ling, and Morawska (2014) found that particles over 10 nm sampled near busy roads were more highly charged than usual. They also showed that the mean concentration of charged particles decreased following a -0.6-th power law of the distance from the roads. Further,  showed that the net charge of PM2.5 collected in a residential area tends to be negatively charged.
While the net charge and concentrations of charged particles have been well studied, there are few studies on the charging states of individual particles in the atmosphere, which are directly related to environmental problems. Very recently,  developed a parallel-plate particle separator (PPPS) that can classify particles by their charging states similarly to the parallel-plate differential mobility classifiers (e.g., Alsharifi and Chen 2018), and measured the charging states of individual particles in ambient aerosol.
In the experiment, however, the PPPS was found to suffer from an inaccuracy in classification: particles flew out to unexpected channels or they remained inside the device. These problems were observed even when a voltage was not applied, where all particles are expected to flow out to the central exit channel. Therefore, we need to elucidate the reason for these problems and to improve the device performance for more accurate measurement.
The objective of the present study is to clarify the mechanism of classification inaccuracy in the PPPS and to provide a solution for performance improvement, by numerical simulation. Another objective is to establish a mathematical method to determine the charge distribution from the number of particles counted at the exit of the PPPS. The rest of this article is organized as follows. The PPPS and its operation are briefly introduced in Section 2. The numerical procedure is outlined in Section 3, followed by the simulation results presented in Section 4. First, the behavior of non-charged particles under the base condition, that is, the same condition as that used in the reference experiment, is investigated. Subsequently, we look for the condition that improves the classification accuracy by altering the flow conditions and the applied voltage. On the basis of the results, a method to determine the charge distribution of aerosol particles using the PPPS is proposed and verified in Section 5; and finally, conclusions are drawn in Section 6. Details of the experimental procedure, verification of the simulation, details of the results of the parametric study, discussion on the differences between the simulation and the experiment, and a discussion on the accuracy of the charge distribution determined by the proposed method are documented in online supplemental information (SI).

Principle of operation of PPPS
The PPPS developed by , called K-MACS (Keio-Measurement system of Aerosol Charging State), consists of a pair of parallel plates (450 mm long, 100 mm high) arranged so that the width between the two plates is 5 mm. Figure 1 shows the schematic of the PPPS at the center of the gravitational direction (i.e., z Ã ¼ 50 mm from the bottom edge). In the figure, all lengths are non-dimensionalized by the half-width between the two plates; that is, h Ã ¼ 2:5 mm. Hereafter, quantities without superscripts denote dimensionless quantities, while dimensional quantities are explicitly denoted by asterisks: this is because the rest of the present analysis is mostly performed in dimensionless form. The PPPS consists of three sections: an inlet section, a classification section, and an exit section.
The inlet and exit sections have three 1-mm-wide flow channels separated by 1-mm-thick parallel plates (see Figure S1 in SI). An air flow of 1.5 L/min is introduced by three external pumps with a flow rate of 0. 5 L/min installed in the downstream of each exit channel. Aerosol (i.e., air with particulate matter) is introduced into the central channel of the inlet section, and sheath air is introduced into the other two inlet channels (see also Figure S2). In the central region of the classification section, electrodes are installed on both sides of the outer plates, and a uniform electric field is imposed between them. As particles flow through the classification section, they gradually drift toward the electrodes by Coulomb's force. Subsequently, each particle flows out to one of the three channels of the exit section depending on the applied voltage and the charging state. Consequently, particles are classified by their charging states.
The number of classified particles is counted by using optical particle counters installed at the downstream of each exit channel. Other properties of classified particles (e.g., chemical composition) can also be measured by installing the appropriate equipment. Additional details on the operation of the PPPS are described in  and SI.

Numerical procedure
The computational domain and the names of different parts are shown in Figure 1. The origin of the xcoordinate is set at the left end of the computational domain, and that of the y-coordinate is at the center of the computational domain. The present simulation is carried out in a two-dimensional field, considering only the central region of the gravitational (z) direction, by assuming that the PPPS is locally uniform in the z-direction. Because the height of the PPPS is 20 times longer than the width between the two electrodes, three-dimensionality such as the edge effect of the flow field can be neglected. Note that Kulkarni and Wang (2006) applied a similar assumption to an analysis of the transfer of aerosol particles in a device formed by two parallel-plate electrodes.

Flow field
The PPPS was designed so that the flow field remains laminar. The governing equations of the flow field are the incompressible continuity equation and Navier-Stokes equation in Cartesian coordinates, expressed as r Á u f ¼ 0; (1) Here, u f is the fluid velocity, p is the pressure, and Re h is the Reynolds number, defined as where U Ã b is the bulk-mean velocity in the classification section, h Ã is the half-width of the classification section, and Ã is the kinematic viscosity of air. The Reynolds number is set at Re h ¼ 7:5 following the experimental condition described in Section 2. Note that all variables without a superscript are non-dimensionalized by U Ã b ; h Ã , and the fluid density q Ã f . The boundary conditions are summarized in Table 1. The labels in the table correspond to those in Figure 1. The inflow and outflow, that is, (a) to (c) in Table 1, are regarded as fully developed laminar flows. Here, the subscript of the bulk-mean velocity denotes the inlet and exit channels. The total flow rate is kept constant among all simulations, although we vary the flow rates of the three inlet channels as parameters in Section 4.2. The boundary condition on the walls ((d) to (f) in Table 1) is the no-slip condition.
The present flow simulation code has been developed on the basis of the direct numerical simulation code for turbulent channel flows (Fukagata, Kasagi, and Koumoutsakos 2006). The governing equations and boundary conditions are discretized over the computational domain on a staggered mesh using the QUICK scheme (Leonard 1979) around the separators at 0 x 10; 153:84 x 163:84 and the secondorder central difference scheme in the other region, that is, 10 x 153:84. Time integration of the unsteady computation is performed by using the third-order Runge-Kutta/Crank-Nicolson scheme with the higher-order SMAC-like velocity-pressure coupling scheme (Dukowicz and Dvinsky 1992) until the flow field becomes steady. After that, particle tracking is performed on that steady flow field.

Electric field
We assume that the charge of the particles does not affect the electric field, because the targeted aerosol is dilute. Therefore, the electric field E is steady and is computed as follows: where U is the electric potential computed by the Laplace equation: These governing equations are non-dimensionalized using the applied voltage / Ã and h Ã . Table 1 shows the boundary conditions. Because the voltage is applied to the electrodes installed on the walls of the center of the classification section, the boundary condition is given as follows: U ¼ 1 at f ð Þ : 40:96 x 122:88 ; y ¼ À1:00; 0 at a ð Þ to e ð Þ in Figure 1: The governing equations and boundary condition are discretized on the same staggered grid as that for the flow field. They are numerically solved using the Gauss-Seidel method until the total sum of residual becomes less than 10 À14 .

Lagrangian particle tracking
Particle trajectories are computed by Lagrangian particle tracking. The equations of motion of a charged particle including Brownian diffusion and the electrostatic force (Ounis, Ahmadi, and McLaughlin 1991) are given by where x p and u p denote the position and velocity of the particle, s p is the particle relaxation time, G is a pair of Gaussian random numbers with zero mean and unit variance, Sc is the Schmidt number, Dt is the computational time step, q is the charge of the particle, and m is the mass of the particle. Again, these equations of motion are non-dimensionalized by U Ã b ; h Ã , q Ã f , and / Ã . The first term on the right-hand side of Equation (8) is the Stokes drag. The particle relaxation time, s p , is given by where d p is the diameter of particle, S is the relative density ratio to fluid, and C C , the Cunningham slip correction factor, which corrects the Stokes drag of submicron particles, is given as Poiseuille flow inlet (Bulk-mean velocity: Poiseuille flow outlet (Bulk-mean velocity: where k Ã is the molecular mean free path of air. The second term on the right-hand side of Equation (8) is the Brownian force. The Schmidt number, Sc, is given by where k Ã is the Boltzmann constant, and T Ã is the temperature. Because the timescale of the Brownian motion is much shorter than the computational time step, the Brownian force is regarded as an independent white noise process in each direction. Note that the amplitude of the Brownian force, i.e., the variance of the white noise, is the same as that described by Fukagata, Zahari, and Bark (2004), who discussed the dynamics of Brownian particles in a turbulent channel flow, and Ramechecandane et al. (2011), who simulated nanoparticle dispersion in DMAs. The third term on the right-hand side of Equation (8) is the Coulomb force. The charge, q Ã , is given as where N c is the number of charges, and Ã is the elementary charge. The mass of a particle, m Ã , is given by assuming spherical particles. Equation (8) is based on the particle equation of motion presented by Maxey and Riley (1983). Ahmadi (1990a, 1990b) showed that the Basset term, Fax en correction, virtual mass term, and pressure term have little effect on the behavior of submicron-scale rigid particles such as those targeted in the present study. Furthermore, Ounis, Ahmadi, and McLaughlin (1991) showed that the Saffman lift is negligible for Brownian particles. Hence, in the present study, these terms are not taken into account. Because the targeted aerosol is regarded as dilute, the present simulations are performed with one-way coupling: the influence of particles on the flow field and particle-particle interactions are neglected.
The initial distribution of particles is set as a uniform distribution at Inlet 2 ((b) in Figure 1), idealizing that particles enter Inlet 2 from the pipe connected vertically to the channel (see also Figure  S1). Collisions between a particle and the walls of the separators are assumed to be perfectly elastic. Once a particle flows out to the exit channels or deposits on the parallel plates with the electrodes, computation of the particle is completed, and its position is recorded ((c), (e), and (f) in Figure 1).
Random numbers with a uniform distribution between 0 and 1 are generated by a portable random number generator with the linear congruential method, called RAN1 (Press et al. 1986), whose period is considered infinite in the practical sense with no sensible sequential correlations. The Gaussian random numbers for G are generated by transforming the uniformly distributed random numbers by using the Box-Muller method (see Press et al. 1986). The velocity and electric fields at the particle locations (u f , E) are obtained by using bilinear interpolation. Time integration is performed by the first-order Euler implicit scheme.
An alternative method for numerical integration of the equations of motion of Brownian particles was developed by Ermak and Buckholz (1980). Their method is based on the analytic solution of the stochastic Langevin equation of a single particle for a small time step. On the other hand, the present method directly simulates the Brownian force as a white noise process that is added to the equation of motion of the Brownian particles. As a result, it is somewhat easier than the method of Ermak and Buckholz (1980) to use a time step that is much larger than the time scale of the Brownian motion and that is suitable for numerical integration of the flow field.

Computational conditions and verifications
The number of grid points for computing the flow field and electric field in the x-direction and the y-direction are 2048 and 200, respectively. The grid resolution has been proved adequate by confirming that the flow field computed with a finer grid, which is 4096 Â 400, has shown acceptable convergence: the maximum difference in the wall-normal velocity v f between the finer grid and the original grid is less than 0.39% at x ¼ 10.24, where both the streamwise velocity u f and the wall-normal velocity v f change most rapidly in the x-direction.
The computational time step Dt is set at 1 Â 10 À3 . The dependence of time steps on particle trajectories and the model of Brownian diffusion have been validated by comparing the variance of y p in a uniform flow with the results obtained by using a shorter time step, Dt ¼ 1 Â 10 À4 , and the exact solution (see Figure S3 in SI). It has been confirmed that the time step is sufficiently small and that the behavior of Brownian particles is in good agreement with the exact solution.

Results and discussion
The air is assumed to be at a temperature of T Ã ¼ 287 K with a density of q Ã f ¼ 1:225 kg=m 3 and a kinematic viscosity of Ã ¼ 1:50 Â 10 À5 m 2 =s. The particle is assumed to be a standard polystyrene particle with a diameter of d Ã p ¼ 0:37 lm and a relative density ratio to fluid of S ¼ q Ã p =q Ã f ¼ 804. The specifications of the particles are summarized in Table 2.

Flow field and particle trajectories in the base condition
In this section, the simulation results under the base condition (hereafter referred to as Case 100) are presented. The widths of the inlet and exit channels and separators are all the same: The bulk-mean velocity in the inlet and exit channels is It is confirmed that the obtained flow field is steady and symmetric with respect to the x-axis after substantial time integration (see Figure S4 in SI). The trajectories of particles are approximately aligned with the streamlines because the particle velocity immediately relaxes to the fluid velocity due to the very short particle relaxation time. However, some particles are observed to flow out to Exits 1 and 3 despite the fact that all the streamlines from Inlet 2 are connected to Exit 2.
In order to clarify the influence of Brownian diffusion, we compare the results in Case 100 computed with and without Brownian diffusion. The probability density functions (PDFs) of the position of particles, p(y), at the end of the computational domain (x ¼ 163.84) are shown in Figure 2. Particles flowing out to Exits 1 and 3 are observed only in the case with Brownian diffusion. Without Brownian diffusion, all particles perfectly follow the streamlines and flow out to Exit 2. This comparison confirms that the outflow of particles to Exits 1 and 3 is due to Brownian diffusion. The PDF has its local maxima at y ' 60:2 in the case without Brownian diffusion. It is found that particles very close to the walls of Inlet 2 flow out from locations slightly away from the walls of Exit 2, where the local maxima are observed. This is because the flow is slightly separated at the corners of the separators due to a sudden contraction. The local maxima disappear in the case with Brownian diffusion because such particles flow out to Exits 1 and 3.
For a quantitative discussion on the device performance, we define the collection ratio C R as where N E;i is the number of particles flowing out to Exit i (i ¼ 1, 2, 3) or that depositing on the walls inside the device (i ¼ 4). Note that the definition of C R is similar to that of the transfer function (Knutson and Whitby 1975), in particular, that defined by Hagwood, Sivathanu, and Mulholland (1999): C R corresponds to the transfer function when it is plotted as a function of the electrical mobility. Table 3 shows C R in Case 100 with and without Brownian diffusion. For Exit 2, to which all particles are expected to flow out in the absence of an electric field, C R;2 ¼ 93:1% with Brownian diffusion, while C R;2 ¼ 100% without Brownian diffusion. Table 3 also shows the mean and standard deviation of C R from 11 data sets obtained using different random seeds. Because the result with the single data set is within the standard deviation of 11 data sets, it is confirmed that a single trial (i.e., tracking of 500; 000 particles) is sufficient for the parametric study conducted below. The results presented below are those obtained by a single trial with the same random seed.

Parametric study for improvement of classification accuracy
In order to suppress the outflow of particles to Exits 1 and 3, we vary the flow rate of Inlet 2 while  maintaining the total flow rate constant and compare the result with experimental observation. The experimental procedure, setup, and equipment are described in SI (see Figures S1 and S2). The cases are referred to as Case Uxxx, where "xxx" represents the flow ratio, F R , defined as the ratio of the flow rate of Inlet 2 to the average of the flow rate of the three inlet channels; that is, The other parameters (e.g., W IC;2 ¼ W IS ¼ 2=5) are unchanged from Case 100. Note that the flow rates of the three exit channels are the same and remain unchanged in all the cases; that is, Figure 3 shows relationship between C R and F R . A dramatic improvement in C R;2 with decreasing F R is observed, and this tendency broadly agrees with the experimental results: in particular, C R;2 is greater than 99% in Case U078 and Case U084. On the other hand, a slight difference in C R is observed between the simulation and the experiment: C R;2 in the simulation is systematically higher than that in the experiment (it is more pronounced for smaller F R;2 ), and C R;3 >C R;1 in the experiment. Table 4 shows C R;2 and the difference in its values for numerical simulation and experiments. The maximum of the difference is 23.23% at the flow ratio F R ¼ 96:2%. This difference is likely to be due to the nonuniform particle distribution at the inlet in the experiment (detailed in Table  S1 and Figure S13 in SI).
Toward development of a better device, the influence of the widths of the inlet channels is studied while maintaining the flow rate in each inlet channel. Case CWxxx indicates that the channel width of Inlet 2, W IC;2 , is changed to xxx% of Case 100 by shifting the positions of the separators, and Case SWxxx indicates that the thickness of the separators, W IS;2 is changed to xxx% of Case 100. It is found that C R;2 increases as W IC;2 decreases and that the C R;2 values of these two cases are similar to each other (detailed in Figures S5 in SI). An approximately 50% improvement is achieved in Case CW025 and in Case SW175.
In order to find the mechanism of improvement, we focus on the particles flowing close to the upper wall of Inlet 2, which tend to flow out to Exit 3. With decreasing F R or W IC;2 , the particle trajectories shift inward following the contraction of streamlines (see Figures S6-S8). Consequently, particle trajectories are more focused toward the desirable exit.
We have also investigated influence of the widths of the exit channels. However, the results show that the change in C R;2 is within twice the standard deviation of the Brownian diffusion shown in Table 3. That is, changing the channel widths in the exit section has virtually no influence on the device performance. This is primarily because the flow rates in the three exit channels are always kept the same in the present configuration following the experimental condition.

Influence of applied voltage on classification accuracy
In this section, we investigate the influence of the width and flow rate of the inlet channel and the applied voltage on the classification accuracy for charged particles. The electric field between the electrodes is nearly uniform in the x-direction and negligibly small elsewhere (see Figure S9 in SI).
The relationship between C R and / Ã is investigated for particles with a constant number of charges, N c ¼ 4 (see Figure S10). This simplification is justified because we have confirmed in our simulation that the movement of an individual particle is identical (except for Brownian diffusion) as long as the product of N c and / Ã (i.e., the terminal velocity in the y-direction calculated by the product of the electrical mobility (Hinds 1982) and the electric field) is the same. In regard to the collection ratio of Exit 3, C R;3 in both Case U084 and Case CW025 is found to be substantially higher than C R;3 in Case 100 (see Figures S10 and S11). The mechanism of the improvement, that  is, the contraction of particles at the entrance of the classification section, is maintained in the presence of an electric field (see Figure S12). In particular, C R;3 in Case U084 is the highest, and C R;3 >99% with / Ã ¼ 130 -135 V. Therefore, the optimal range of the applied voltage to classify particles with N c ¼ 4 is considered to be / Ã ¼ 130 -135 V. Because it has been confirmed that the behavior of charged particles can be scaled by the electrical mobility, the behavior of particles charged with an arbitrary number of charges N c (hereafter referred to as N ccharged particles) can be extrapolated by using the results of N c ¼ 4. The optimal voltage / Ã N c for N ccharged particles can therefore be expressed as 520 Figure 4 shows the C R;3 of N c -charged particles computed using the results of N c ¼ 4. It is shown that the negatively charged particles never enter Exit 3 under a positive voltage. However, we can see that it is difficult to collect the particles with a particular charge number only, because the lines of C R;3 for different N c values overlap with each other.
For classification under different conditions (e.g., width between the two plates, particle diameter, and flow rate), non-dimensional analysis of the classification is conducted. Regardless of the conditions, the traveling distance in the y-direction non-dimensionalized by h Ã , that is, y e ¼ y Ã e =h Ã , is considered to be constant. Here, y Ã e is calculated by the electrical mobility, Z Ã p , and the residence time of particles inside the electrodes, t Ã e , as where Z Ã p is defined as Because t Ã e can be estimated from the bulk-mean velocity U Ã b and the length of the electrodes, L Ã , as Finally, y e ¼ y Ã e =h Ã can be rearranged as Here, Z Ã p E Ã =U Ã b represents the relative velocity toward the electrodes, and L Ã =h Ã is the relative length of the electrodes. Figure 5a shows C R;3 as a function of y e in Case U084. The figure corresponds to a transfer function of the PPPS for general conditions. Here, it is found that C R;3 takes its maximum of C R;3 ¼ 100% at y e ¼ 0:72. Moreover, Figure 5b shows the edges of C R;3 (i.e., C R;3 ¼ 0%; 100%) under various non-dimensional conditions decomposed on the vertical axis of Z Ã p E Ã =U Ã b and the horizontal axis of L Ã =h Ã . These results provide a prediction for C R;3 and the appropriate applied voltage when the experimental conditions and configuration of electrodes of a PPPS-like device are substituted.

A proposal to measure the charge distribution of aerosols using PPPS
On the basis of the above results, we propose and demonstrate a method for measuring the charge distribution of aerosols from the numbers of particles counted at the exit section of the PPPS.

Proposed method
We consider an aerosol that consists of n i particles with a charge number of N c ¼ i. Under a voltage / Ã , the number of particles counted at Exit 1, N E;1 ð/ Ã Þ, can be expressed as where C i R;1 ð/ Ã Þ is the collection ratio of particles with  Thus, one can determine n i by solving a system of linear equations of rank equal to the number of unknowns. Because N E;1 ð/ Ã Þ ¼ N E;3 ðÀ/ Ã Þ due to the symmetricity, the actual number of applied voltages required for determining the unknowns is halved. The equations at all exits and applied voltages are represented in a matrix-vector form as Here, the matrix C R can be obtained by reading Figure 4. Therefore, one can simply determine n by multiplying the inverse matrix of C R as Finally, the PDF of N c -charged particles, pðN c Þ (i.e., the charge distribution) can be computed by normalizing n by the total number. It should be noted that a set of / Ã must be chosen so as to retain the rank of C R .

Validation
As a test, we attempt to determine the charge distribution of particles using this method. The Boltzmann charge distribution (Cooper and Reist 1973), which is the common charge distribution of aerosol in the air, with a range of À5 N c þ5, is used as the charge distribution of particles introduced from Inlet 2. The applied voltages are selected as 108; 135; 180; 270; and 540 V, which correspond to the optimal voltages for the different charge numbers (Equation (16)). Figure 6 shows the charge distribution of particles, pðN c Þ, determined using the proposed method by counting the numbers of particles at exits under different voltages, N E . The computed charge distribution is found to be in excellent agreement with the input charge distribution. We define the residuals R i as where the superscript ðÁÞ i denotes the charge number. It is found that R i reaches a maximum of R 0 ¼ 5:2 Â 10 À3 , and the sum of R i is 1:5 Â 10 À2 . These results demonstrate that the proposed method is valid for determining the charge distribution of aerosol particles using the PPPS. The results of the tests under various conditions (summarized in Table S2) are discussed in SI. It is shown that the accuracy of the determined charge distribution depends on the determinant of C R , which is determined by the flow field and applied voltages, and therefore selection of an appropriate set of applied voltages is also important for the estimation accuracy in addition to the improvement of flow as discussed in Section 3 (see Figure S14). In addition, it is shown that the proposed method is reproducible in experiments because the charge distribution determined with significant figures of two digits indicates allowable errors (see Figures S15 and S16).
It is also worth noting for its practical use that the present method can be extended to the case where the number of Equation (22), that is, the number of measurement times, is larger than the number of elements of n, that is, the number of unknowns, by solving the Moore-Penrose pseudoinverse matrix of C R (e.g., pinv in MATLABV R ) instead of C À1 R . Further, solvers for the non-negative linear least-square problem (e.g., lsqnonneg in MATLAB V R ) are useful when the probability density pðN c Þ obtained by Equation (24) becomes negative due to errors in measurement.

Conclusions
The trajectories of submicron particles in a PPPS are numerically analyzed using the Lagrangian particle tracking method taking into account the Brownian force and the electrostatic force in order to clarify the mechanism of the classification uncertainty and to improve the device performance.
First, it has been confirmed that Brownian diffusion decreases the classification accuracy of the PPPS: 7-8% particles flow out to unexpected exits due to Brownian diffusion even in the absence of an electric field. The device performance is improved by decreasing the flow rate of the aerosol inlet, which is in accordance with the experimental observations. Changing the widths of inlet channels, by both narrowing the aerosol inlet and widening the separators, is also found to substantially improve the device performance. The improvement is caused by the contraction of flow, on which particles trace before they enter the classification section. The above results provide a useful guideline for a better design or flow condition of the PPPS. Further, a non-dimensional analysis of the classification of particles is conducted. The result enables the prediction of classification (i.e., transfer function) under various experimental conditions and configurations of electrodes.
Finally, on the basis of the simulation results, we have proposed a method to determine a charge distribution of aerosol from the numbers of particles counted at the exit section of the PPPS. It has been demonstrated that the proposed method can perfectly recover the charge distribution given at the aerosol inlet from the number count at the exits.

Nomenclature
Latin symbols C C Cunningham slip correction factor [N.D.] C R collection ratio defined as Equation (14)  Superscripts and subscripts (no superscript) non-dimensional quantity ðÁÞ Ã dimensional quantity Inlet i or Exit i (for i ¼ 1; 2; 3) or walls inside the PPPS (for i ¼ 4) ðÁÞ p quantity for particles Abbreviations DMA differential mobility analyzer K-MACS Keio-Measurement system of Aerosol Charging State PDF probability density function PPPS parallel-plate particle separator