The Effect of Nanoparticle Convection-Diffusion Loss on the Transfer Function of an Aerosol Particle Mass Analyzer

Abstract

The existing theoretical response spectra of APM-3600 agree well with the experimental data for submicron particles larger than 100 nm in the electrical mobility diameter but not for nanoparticles. In this study, a 2-D numerical model was developed to predict the transfer function and response spectra of APM-3600 based on the detailed simulation of flow and particle concentration fields. It was found that recirculation flows existed in the annular classifying region and APM's inlet and outlet regions, which led to enhanced convection-diffusion loss of nanoparticles compared to that without considering flow recirculation. As a result, the APM underestimates the mass of naonoparticles due to the shift of the peak position of the transfer function to a larger diameter than the targeted diameter. The response spectra calculated with the simulated transfer function agree well both in shapes and peak values with the experimental data present in a previous study for both nanoparticles and submicron particles larger than 100 nm. The predicted particle masses also agree well with the PSL's experimental data of the article.

Copyright 2014 American Association for Aerosol Research

INTRODUCTION

There have been extensive laboratory and ambient air studies related to particle mass, mass concentration, density, and morphology (McMurry et al. 2002; Park et al. 2003a, 2004; DeCarlo et al. 2004; Geller et al. 2006; Virtanen et al. 2006; Lall et al. 2008; Lee et al. 2009; Tsai et al. 2009, 2011; Chen et al. 2010, 2013; Awasthi et al. 2013; Johnson et al. 2013; Symonds et al. 2013). Two commercially available particle mass analyzers, namely, aerosol particle mass analyzer (APM) (Ehara et al. 1995, 1996) and Couette centrifugal particle mass analyzer (CPMA) (Olfert and Collings 2005; Olfert et al. 2006) were often used in these studies. Both classify particles based on the mass-to-charge ratio or specific mass, S (kg/C), which is defined as [1] where m and q are the mass (kg) and charge (C) of the particle, respectively.

The APM classifying region is the space between two coaxial cylinders in which charged particles are classified by the centrifugal force (F_c) created by the cylinders rotating at the same angular velocity and the electrostatic force (F_E) in the opposite direction created by a voltage drop between the cylinders (Ehara et al. 1995, 1996). The CPMA operates with the same principle as that of the APM except that the inner cylinder rotates slightly faster than the outer cylinder to improve the particle transmission efficiency (Olfert and Collings 2005). In this study, only the transfer function and the response spectra of APM-3600 are investigated, but the results of the study would also work for the CPMA or other model of the APM.

The transfer function of particles with a specific mass S, Ω_APM(S), is the probability for the particles to pass through the classifying region (Ehara et al. 1995, 1996; Hagwood et al. 1995; Olfert and Collings 2005) or the whole APM (in the study). Therefore, the transfer function defines the performance of the APM in which the peak position of the transfer function determines the accuracy of the mass measurement, the width determines the resolution, and the height determines the particle loss in the APM. Hence, an accurate model of the transfer function is crucial to the accuracy of the APM. Currently, the model developed by Ehara et al. (1995, 1996), or called Ehara model, is the most commonly used model in which the particle mass is calculated as [2] where m_c is the central particle mass or the targeted particle mass (kg), which was considered as the measured mass; q is the number of elementary charges; e is the elementary charge (1.6 × 10⁻¹⁹ C); V and ω are the voltage and angular rotational speed of the APM (rad/s), respectively; r₁ and r₂ are the inner and outer radii of the classifying region (m), respectively; r_c is the average of r₁ and r₂ (m). Ehara model is accurate for submicron particles. However, it was reported that both APM-3600 and APM-3601 (the compacted APM) underestimated the mass of monodisperse nanoparticles smaller than 30 nm (PSL, NaCl, Santovac) (Tajima et al. 2011, 2013). Brownian motion and morphology of particles in the APM were suspected as the causes of underestimation.

The response spectra of the APM, the ratios of total particle number concentration at the APM outlet, N_OUT (#/m³), to that at the APM inlet, N_IN (#/cm³), can be calculated as the products of the particle-size distribution at the inlet of the classifying region and the transfer function as [3] where S_min and S_max are the minimum and maximum specific masses that can pass through the APM; N_in(S) is the mass distribution of the particles at the APM inlet. Since the transfer function cannot be verified by the experimental data directly because of perfectly monodisperse particles cannot be generated experimentally, the response spectra are usually used to verify the models of the transfer function. Some studies assumed N_in(S) to be a δ function (Ehara et al. 1995), perfectly monodisperse distribution (Olfert et al. 2006), narrow triangle size distribution (Lall et al. 2009), and Gaussian distribution (Tajima et al. 2011) when calculating the response spectra.

Several studies have developed models to predict the transfer functions as well as the response spectra of the APM (Ehara et al. 1995, 1996; Hagwood et al. 1995; Olfert and Collings 2005). Again, Ehara model is valid for submicron particles but overestimates the response spectra of nanoparticles smaller than 50 nm (Ehara et al. 1996; Tajima et al. 2011, 2013). Hagwood et al. (1995) presented a stochastic differential equations (SDE approach) and a Monte Carlo model (MD approach), respectively, to predict the transfer function of a prototype APM taking Brownian motion of particles into account. Good agreement of predicted transfer functions between two approaches was obtained. For the rotational speed of 3000 rpm and the voltage of 0.26 volt of the APM, the peak value of the theoretical transfer function of particles of 20 nm in diameter was found to decrease from about 0.9 to 0.2 when the diffusion was taken into consideration in the SDE approach. It was also found that the peak of the transfer function was shifted to a larger diameter than the targeted diameter of 20 nm but the reason was not clearly explained.

In addition, the MD approach was found to over-predict the peak height of the transfer function of the APM (APM-10, Kanomax Inc., Japan) for particles with d_p = 300 nm due to diffusion and transport losses inside the APM and during the transport to the CPC (Lall et al. 2009). Therefore, the applicability of the SDE and MD approach needs to be further examined for nanoparticles. Olfert and Collings (2005) developed an Eulerian diffusion model with convection-diffusion equations to numerically calculate the transfer function of the APM and CPMA. The amplitude of the experimental transfer function was found to be lower than the simulated value, which might be caused by the particles lost in the analyzer in the regions before and after the classifying region (Olfert et al. 2006).

The objective of this study is to develop numerical models to predict the transfer function of nanoparticles accurately with detailed simulation of flow and particle concentration fields. The models will also be validated with experimental data for particle mass and response spectra in Tajima et al. (2011).

NUMERICAL METHOD

Two numerical models, model 1 and model 2, were used to predict the transfer function of APM-3600 at the fixed aerosol flow rate of 1 L/min. Model 1 is based on the parabolic flow profile, in which only the transfer function in the classifying region was predicted as in previous studies (Ehara et al. 1995, 1996; Hagwood et al. 1995; Olfert and Collings 2005). Model 2 is based on the detailed flow profile, in which the transfer function in the whole region (inlet, outlet, and classifying regions) was predicted. Calculation domain is shown in Figure 1 in which the regions covered by thick lines are rotating. r₁ and r₂ are the inner radius (0.05 m) and outer radius (0.052 m) of the classifying region, respectively; r₃ is the radius of the inlet and outlet tubes of APM-3600; L is the length of the classifying region (0.25 m); l₁, l₂, l₃, l₄, l₅, and l₆ are regions before and after the classifying region. In Figure 1, different scales of 1:12.6, 1:18.3, 1:240, 1:6.88 are used for the r₁, r₂, r₃ and L, while 1:5.2, 1:3.25, 1:96.7, 1:260, 1:1.45, and 1:1.27 are used for regions l₁ to l₆, respectively. Otherwise some regions will become invisible.

Figure 1 The schematic diagram of APM-3600. The regions covered by thick solid lines are rotating.

A total of 7564 (122 in the direction or the axial direction, 62 in the r direction or the radial direction) and 13,104 (182 in the direction, 72 in the r direction) nonuniform rectangular girds were used in models 1 and 2, respectively. The average grid size is about 2 mm and 32.3 μm in the and r direction, respectively, with the smallest size of 3.6 μm near the wall of the cylinders in model 1; and 5.4 mm and 722 μm in the and r direction, respectively, with the smallest size of 125 μm near the wall of the cylinders in model 2. It was found that increasing the number of girds from 7564 (122 × 62) to 10,004 (122 × 82), and 13,104 (182 × 72) to 16,744 (182 × 92) in models 1 and 2, respectively, only resulted in a slight increase in predicted particle diffusion loss. For example, the diffusion loss of 30 nm particles only slightly increases from 82.5% to 82.6% and 69.5% to 69.6%, respectively. Thus, fixed numbers of grids of 7564 and 13,104 were used for model 1 and 2, respectively, with the total number of iterations of about 100,000 and 4,000,000, respectively, to reach the convergence for the flow field.

Flow Field

For model 1, the parabolic flow field in the classifying region is given by [4] [5] where and _max are the flow and the maximum flow velocities in the direction (m/s), respectively; Q is the aerosol flow rate (m³/s); r is the radial position of the particle (m).

For model 2, laminar flow is assumed because the maximum axial Reynolds number (Re_max) based on the diameter of the inlet tube or outlet tube is about 685, which is much smaller than 2000 (Hinds 1999). The 2-D Navier–Stokes equations are [6] [7] where u and are the air velocities in the r and direction (m/s), respectively; is the air velocity in the θ direction, which is equal to rω; ρ_a is the air density (kg/m³); P is the pressure (Pa); μ_a is the air dynamic viscosity (kg/m·s). The flow is assumed to be incompressible and the continuity equation is [8]

The Navier–Stokes and continuity equations were discretized by the finite volume method and solved by the SIMPLER algorithm (semi-implicit method for pressure-linked equations) (Patankar 1980) with the following initial setting of velocities: [9] [10] [11] [12]

The boundary conditions are: [13] [14]

Particle Concentration Field

The governing equation for the particle concentration, N_p, is [15] where τ = ρ_pd_pC_c/18 μ_a, which is the particle relaxation time (s); ρ_p is the particle density (kg/m³); C_c is the Cunningham correction coefficient (Hinds 1999); Z_p is the particle electrical mobility (m²/s·V), which is defined as Z_p = qeC_c/3πμ_ad_p; D_B is the Brownian diffusion coefficient for particles (m²/s); E_r is the electric field strength (volt/m) in the classifying region, which is given by [16]

Perfect absorption was assumed at the walls where particle concentrations were zero. Equation (15)[15] was also discretized by using the finite volume method and solved by the same computer code as that used in the flow field simulation. About 10,000–30,000 iterations were needed to reach convergences. After the profiles of the flow field and particle concentration were obtained, the transfer functions of model 1 and model 2, which are Ω₁ and Ω₂, respectively, were calculated as [17] [18] where L_total is the total length of the APM, or the sum of l₁, l₂, l₃, L, l₄, l_5, and l₆. Particles were assumed to be spherical with a known density and singly charged. Therefore, the mass-to-charge ratio can be converted to particle diameter readily.

RESULTS AND DISCUSSION

Particle Penetration of Still APM-3600

Particle loss in still APM-3600 (0 rpm, 0 volt) is mainly due to diffusion loss (Tajima et al. 2011). Model 1, model 2, and Gormley and Kennedy equation (Hinds 1999) were applied to predict the penetration of particles. Particle penetration based on Gormley and Kennedy equation, P_G&K, can be written as: [19a] [19b] where the dimensionless deposition parameter μ is [19c]

Predicted particle penetrations are compared with the experimental data of Tajima et al. (2011) as shown in Figure 2. The predicted penetrations by model 1 and P_G&K agree very well with a deviation less than 0.6% for d_p less than 200 nm. The agreement between model 1 predictions and the experimental data is good for d_p > 50 nm, but overestimation exists as much as 14%–26.5% for nanoparticles less than 30 nm. It is because only the particle penetration in the classifying region is predicted in both model 1 and P_G&K, the particle diffusion loss in regions l₁–l₃ and l₄–l₆ is neglected. In contrast, model 2 considers all regions of the APM, and the predictions agree well with the experimental data with a deviation less than 10% for 10 nm ≤ d_p ≤ 200 nm.

Figure 2 The experimental and predicted particle penetrations in still APM-3600 with zero rotational speed and zero voltage.

To elucidate the particle diffusion loss in regions l₁–l₃ and l₄–l₆, the calculated particle concentration contours of 30-nm particles at zero rotational speed and zero applied voltage is shown in Figure 3a as an example. It is noted that a mismatch of contours at the upper right hand side near the classifying region and region l₄ is due to different scales used in different APM regions. Particle concentration contours in regions l₃ and l₄ are too small to visualize if the same scales as that in the classifying region are used. The normalized particle concentration close to the boundaries of the APM seem to be reduced near the walls of regions l₁–l₃ and l₄–l₆ due to convection-diffusion loss. In addition, the recirculation flow is observed in region l₃, as shown in Figure 3b, which enhances the convection-diffusion loss.

Figure 3 (a) Particle concentration contours and (b) flow streamlines in the still APM. A mismatch of contours at the upper right hand side of the classifying region is due to different scales used in the plot.

Good agreement between the predicted and experimental particle penetrations in the still APM confirms that considerable convection-diffusion loss of nanoaprticle occurs in the inlet and outlet regions. The accuracy of the predicted convection-diffusion loss in the APM is significantly improved by model 2, which considers all regions of the APM with detailed flow field simulation.

Comparison of Transfer Functions

The numerical transfer functions predicted by model 2 are compared with those by model 1 and Ehara model as shown in Figure 4. In the simulation, the dimensionless number, λ_c, which determines the height and resolution of the transfer function (Ehara et al. 1995, 1996) was fixed at 0.22. λ_c can be written as [20] where τ_c is the relaxation time of the center particle with the targeted specific mass or so-called center specific mass (Ehara et al. 1995); ω_c is the central rotational speed corresponding to the targeted particle with the diameter d_pc.

In Figure 4, the predicted transfer functions are for particles with d_pc of 208 nm, 100 nm, 51 nm, 30.6 nm, 20 nm, and 10 nm with given densities. The central rotational speed is calculated with Equation (20)[20] based on the fixed λ_c of 0.22, while the central voltage, V_c, can be derived from the balance between the electrostatic force and centrifugal force of singly charged center particles as [21]

Figure 4 The transfer functions of Ehara model, model 1, and model 2 for (a) 208-nm, (b) 100-nm, (c) 51-nm, (d) 30.6-nm, (e) 20-nm, and (f) 10-nm particle (λ_c = 0.22).

In Figure 4, the gray and black lines are the transfer functions predicted by model 1 and model 2, respectively, while the dash lines are those predicted by the analytical Ehara Model. Since Ehara model neglects convection-diffusion loss, it is considered as the ideal transfer function. The vertical thin black line locating at = 1 corresponds to the diameter (or specific mass) of the targeted or center particle. For 208-nm and 100-nm particles, the transfer functions of all models are very close in both the height and width (Figures 4a–b). However, for particles smaller than 100 nm, the transfer functions of model 2 become lower than those of model 1 (Figures 4c–f), while the transfer functions of Ehara model maintain at the same height of 0.8. This comparison shows the validity of the analytical Ehara model for particles greater than 100 nm, but not for particles smaller than 100 nm.

On the other hand, although model 1 considers the convection-diffusion loss in the classifying region, it underestimates the convection-diffusion loss and the degree of the underestimation increases with decreasing particle diameter (Figures 4c–f). For example, the peak height of the transfer function of 30.6-nm particles at = 1 for model 1 is higher than that of model 2 by as much as 65% (Figure 4d). This is because model 1 considers the particle convection-diffusion loss in the classifying region only, while the loss can also occur in other regions as shown in the previous section. The normalized particle concentration contours shown in Figure 5a for 30.6-nm particles in APM-3600 operating at 6406 rpm and 4.79 volt further indicate the enhanced convection-diffusion loss due to recirculation flow regions formed at such a high rotational speed as shown in Figure 5b. More recirculation regions with stronger recirculation strength occur as compared to those in the still APM (Figure 3b), which lead to enhanced convection-diffusion loss. This causes more nanoparticle loss as predicted by model 2.

Figure 5 (a) Particle concentration distribution of 30.6-nm particles at the applied voltage of 4.55 V. (b) Flow streamlines at the rotational speed of 6406 rpm.

Furthermore, the convection-diffusion loss also increases the width of the transfer function as shown in Figure 4. The calculated transfer functions of both model 1 and model 2 become wider than those of Ehara model for nanoparticles. It is also found that the transfer functions are shifted to the right of the vertical line () due to more convection-diffusion loss of smaller nanoparticles.

The transfer functions of model 1 and model 2 are applied to calculate the response spectra which are then compared with the experimental data of Tajima et al. (2011), who applied APM-3600 to measure the mass distribution of monodisperse PSL particles with a fixed rotational speed but different voltages. The experimental response spectra are the ratios of the measured particle concentrations at the APM outlet to those at the APM inlet. For the simulated response spectra, particles at the APM inlet are assumed to be perfectly monodisperse, which make the response spectra equal to the transfer function. The assumption will not lead to significant errors if the particles used in the experiments are nearly perfectly monodisperse.

Figure 6 shows the comparison of the response spectra for a fixed λ_c value of 0.22. The symbols are experimental response spectra of Tajima et al. (2011), while the thick gray lines and black lines are the response spectra predicted by model 1 and model 2, respectively. To quantify the width of the transfer function, the voltage was normalized by the corresponding V_c, which was calculated with Equation (21)[21] based on the mobility diameter determined by a differential mobility analyzer (DMA) and the bulk densities of particles provided by the manufacturer (JSR Corp., Japan), which were 1.054, 1.06, 1.061, and 1.115 g/cm³, respectively, for 208 nm, 100 nm, 51 nm, and 30.6 nm PSL particles (Tajima et al. 2011). In addition, a vertical thin black line was located at = 1 to indicate the peak position of the response spectra. For 208 nm particles (Figure 6a), the response spectra of both models 1 and 2 agree well with the experimental data. In contrast, for particles equal to or less than 100 nm (Figures 6b–d), the heights of the response spectra of model 2 are closer to the experimental data compared to those of model 1. The absolute differences between the peak heights of the response spectra of model 2 and the experimental values are 6.9%, 4.5%, and 8.3% for 100-nm, 51-nm, and 30.6-nm particles, respectively, which are substantially smaller compared to the larger overestimations of 14.1%, 20.8%, and 22.6% by model 1, respectively. That is, model 2 is more accurate than model 1 for predicting the response spectra of nanoparticles.

Figure 6 The response spectra of (a) 208 nm, (b) 100 nm, (c) 51 nm, and (d) 30.6 nm PSL in APM-3600 operating at λ_c of 0.22.

As shown in Figure 6, the widths of simulated response spectra are seen to be somewhat narrower than those of experimental data especially for smaller nanoparticles. This could be due to uncertainties in the experiment. First, the particles at the APM inlet were not perfectly monodisperse. A particle size distribution with a certain geometric standard deviation would result in somewhat wider response spectra than those predicted here. Moreover, multiply charged particles (e.g., doublets or triplets) with acceptable mass-to-charge ratios can also pass through the APM resulting in wider experimental response spectra (Johnson et al. 2013; Symonds et al. 2013; Tajima et al. 2013). If the experimental size distribution of particles at the APM inlet is used for the prediction of the response spectra, the agreement between the predictions and experimental data will be improved further. In addition, the effect of high rotational speeds (e.g., 6406 rpm for 30.6 nm particles) on the flow field in the θ direction, which was neglected in the study, could also result in the discrepancies due to the potential recirculation flow in the θ direction.

Accuracy of Particle Mass Measurement

An important effect of the enhanced convection-diffusion loss of nanoparticles is on the mass measurement accuracy due to the shift of the transfer function to the right of = 1 as predicted by model 2. For example as shown in Figure 4f, the peak transfer function of 10-nm particles is shifted from 10 nm to 11 nm. A small difference of 9% in the diameter will result in a significant difference of 25% in the particle mass. The shift of the transfer function for 20-nm particles in a prototype APM was also observed by Hagwood et al. (1995) for simulations considering Brownian motion of particles. Tajima et al. (2011, 2013) also mentioned that the APM underestimated the mass of particles smaller than 50 nm. Since the shift of the transfer function predicted by model 1 is not as significant as that by model 2, the enhanced convection-diffusion loss of nanoaprticles is, therefore, the major reason that induces the shift in the transfer function.

The accuracy of the predicted mass, m_pred, can be observed by comparing m_pred/m_c with m_mea/m_c as shown in Figures 7, where m_mea is the measured mass of PSL in Tajima et al. (2011). In the calculation of particle mass, the bulk densities of PSL given by JSR Corp. were used and were assumed to be 1.05 g/cm³ for particles whose densities were not provided_. It is seen that m_pred/m_c agrees with m_mea/m_c indicating model 2 is accurate. For nanoparticles larger than 50 nm, m_pred agrees with m_c or m_pred/m_c approaches 1.0 for both λ_c values of 0.22 and 0.49 as shown in Figures 7a and b, respectively, while it is lower than m_c for smaller nanoparticles (or m_pred/m_c <1.0) and m_pred/m_c ratio decreases with decreasing nanoparticle diameter. The underestimation of measured particle mass is due to the enhanced convection-diffusion loss as mentioned before.

Two different rotational speeds, 3340 rpm and 8465 rpm, were further used to examine the effect of the recirculation flow on transfer functions. The densities of the particles, which were chosen based on the fixed λ_c of 0.22 and the rotational speeds, are listed in Table 1. Figure 8a shows the ratio of the height of the transfer functions (denoted as TFH) for APM-3600 operating at 8465 rpm to that at 3340 rpm. It is observed that the ratio decreases from 1 to 0.76 for a decreasing particle diameter from 100 nm to 10 nm due to stronger recirculation flow (and more nanoparticle loss) at 8465 rpm than that at 3340 rpm. Figure 8b shows the predicted mass over central mass (m_pred/m_c) versus particle diameter at both rotational speeds. The degree of the underestimation is seen to increase with decreasing nanoparticle diameter and increasing rotational speed. The underestimation of mass measurement for 8465 rpm, 25%, is much higher than that for 3340 rpm, 11%, for 10-nm particles. The underestimation becomes less significant for particles larger than 50–60 nm for both rotational speeds with a deviation less than 3%.

Figure 7 The ratio of m_mea to m_act or m_pred to m_act versus the particle diameter. λ_c is fixed at (a) 0.22 and (b) 0.49. The crosses and open circles are the measured number-averaged and mode masses, respectively.

Table 1 Parameters used in Figure 8

APM-3600 (λc = 0.22)
	for 3340 rpm		for 8465 rpm
dp (nm)	ρp (g/cm^3)	m (fg)	ρp (g/cm^3)	m (fg)
100	1.06	0.5550	0.165	0.0864
48	2.51	0.1453	0.39	0.0226
35	3.55	0.0797	0.553	0.0124
30	4.19	0.0592	0.652	0.0092
25	5.08	0.0416	0.79	0.0065
20	6.46	0.0271	1	0.0042
15	8.65	0.0153	1.35	0.0024
10	13.1	0.0069	2.043	0.0011

Figure 8 The ratio of (a) TFH for 8465 rpm to 3340 rpm and (b) m_pred to m_act for 8465 rpm and 3340 rpm versus the particle diameter. λ_c is fixed at 0.22.

Therefore, the enhanced convection-diffusion loss has a large effect on the transfer function of nanoparticles resulting in the underestimation of measured nanoparticle mass when applying Ehara model to determine the nanoparticle mass. Hence, it is recommended that the present model 2 be used and the shift of the transfer functions for nanoparticles less than 50 nm be considered to avoid the underestimation of measured particle mass as reported in Tajima et al. (2011, 2013).

CONCLUSION

This study developed a detailed 2-D mathematical model to predict detailed flow and particle concentration fields in APM-3600. Present model 2 predicts accurately the penetration of nanoparticle in the APM with zero rotational speed and zero voltage. It is shown that diffusion loss of nanoparticles occurs not only in the classifying region but also other regions in the APM. In addition, the convection-diffusion loss reduces the peak and increases the width of the transfer function of nanoparticles. The comparison of the peak heights of the predicted response spectra and the experimental data in Tajima et al. (2011) shows good agreement with a deviation less than 8.0% for particles with 30.6 ≤ d_p ≤ 2 08 nm. The peak position of the transfer function is shifted due to the convection-diffusion loss leading to the underestimation of nanoparticles mass based on Ehara model. In summary, the peak, width, and peak position of the transfer function of nanoparticles are significantly affected by the enhanced convection-diffusion loss by recirculation flows in the APM, which have not been found in previous literatures.

The study concludes that Ehara model is accurate enough for the transfer functions of submicron particles larger than 100 nm, while present model 2 is recommended for more accurate prediction of the transfer functions of nanoparticles. In the future, it is important to further validate the predicted transfer functions considering particle morphology and multiple charge effects for accurate mass and mass distribution measurements of nanoparticles.

FUNDING

The support of Taiwan Institute of Occupational Safety and Health, Council of Labor Affairs (contract number IOSH 102-H324) and National Chiao Tung University (contract number 102W985) is greatly acknowledged.