JOURNAL OF ANALYTICAL ATOMIC SPECTROMETRY JUNE 1994 VOL. 9 701 Simulation of Nebulization Process in Inductively Coupled Plasma Atomic Emission Spectrometry with a Modified Model Using the Monte Carlo Technique* Hu Yanping and Zhang Zhanxiat Department of Chemistry Zhongshan University Guangzhou 5 70275 China Zheng Jianguo Guangdong Commodity Inspection Bureau Guangzhou China An improved Monte Carlo program based on a jet model for the simulation of various particle loss processes that occur in the spray chamber is presented. This program is evaluated by comparing the simulated mass transport efficiency and total mass transport rate with those found by experiment. The aerosol flow velocity distribution as a function of axial and radial distances of the chamber is outlined. The effect of carrier gas flow rate on the distribution of loss particles along the chamber due to various loss processes is studied.The results indicate that the predominant loss mechanism depends strongly on the magnitude of the carrier gas flow rate. Keywords Nebulization process; inductively coupled plasma atomic emission spectrometry; Monte Carlo simulation Generally the liquid samples are converted into aerosols before introduction into the inductively coupled plasma (ICP) through a pneumatic nebulization-spray chamber system. The solution is split into droplets under the influence of a high- speed gas flow and the aerosol produced is highly polydispers- ive with droplet diameters of up to l o o p . Nukiyama and Tanasawa’ have derived empirically the size of the primary aerosol droplets relative to the operating parameters of the nebulizer and the physical properties of the solution.However the equation is valid only for describing the aerosol formed at the nebulizer. Because of other modifying processes,2 namely gravitational settling impaction turbulence and centrifugation taking place in the spray chamber and in the transport of the aerosol between the nebulizer and the plasma the aerosol will have entirely different properties from the primary aerosol on reaching the plasma. These processes generally act to shift the aerosol size distribution to smaller droplets. The spray chamber can be regarded as a hypothetical filter having a cut-off diameter d through which the aerosol is passed i.e. droplets having diameters larger than d are retained and drained to waste while those droplets having diameters of d or less pass through unaffected.Different mechanisms have been suggested to describe the primary separation process in an aerosol chamber. Skogerboe and Olson3 studied gravitational settling and inertial deposition and concluded that gravitational set- tling imposed the primary limitation on aerosol transport for the chamber and condition used. Browner et a1.’ studied gravitational settling impaction turbulence and centrifugal loss and proposed that turbulence-induced loss is the predomi- nant mechanism for the separation of large droplets in an ICP aerosol chamber of the dual concentric type. Gustavsson4 considered that the most likely process causing separation of large droplets in an aerosol chamber is inertial dep~sition.~ Sharp5g6 summarized from the experimental data that for a fixed nebulizer geometry the liquid-to-gas ratio is the param- eter that most determines the quality of the aerosol produced.He further confirmed that the nature of the processes determin- ing droplet deposition rates is particle size dependent and that the inertial and turbulence losses are important whereas the gravitational and centrifugal losses are not important. The aerosol drop size distribution can give important infor- f To whom correspondence should be addressed. * Presented at the XXVIII Colloquium Spectroscopicum Internationale (CSI) York UK June 29-July 4 1993. mation on the sizes of the droplets entering the plasma prior to evaporation decomposition and atomization.The mass transport efficiency defined as the percentage of the mass of nebulized solution that actually reaches the plasma and the total analyte mass transport rate defined as the total mass of analyte reaching the plasma per second are important param- eters in the quantitative characterization of critical aerosol properties as well as in optimizing the experimental param- eters. A further study using a more powerful method is thus essential to understanding all the processes occurring in the spray chamber and to determining the predominant processes responsible for the loss of aerosol droplets. In a previous paper,’ the Monte Carlo technique was used for the first time to simulate nebulization processes in ICP atomic emission spectrometry (ICP-AES).The results were encouraging. The study was based on the fact that the drop size distribution is Gaussian and that the aerosol droplet loss due to gravitational settling impaction turbulence and centri- fugation are random in character. The operating conditions and the effect of the design of the nebulizer on the mass transport efficiency the analyte mass transport rate and the aerosol drop size distribution were obtained in a straight- forward manner since the Monte Carlo technique is especially useful for complex situations and has the ability to obtain information which cannot be extracted experimentally. However the model used previously was too simple. In this work a modified model namely the jet model is proposed with the aim of obtaining more realistic information of the nebulization process and discussing the loss mechanism taking place in the spray chamber.Theory A Monte Carlo simulation is based on the principle that any complex process can be broken down into a series of simpler independent events each represented by a probability distri- bution. Thus the various processes occurring within the spray chamber that generally act to shift the aerosol drop size distribution to smaller droplets can be considered indepen- dently when the time period being viewed is sufficiently small. According to Browner et aZ.? the gravitational turbulence and centrifugal loss processes responsible for the separation of large droplets in the spray chamber can be expressed as follows. (i) Impaction loss dci= 1.5 q g W P ~ V ( 1 )702 JOURNAL OF ANALYTICAL ATOMIC SPECTROMETRY JUNE 1994 VOL.9 where dCi is the cut-off diameter (pm) of the aerosol droplet due to impaction loss q is the absolute gas viscosity ( P ) W is the width of the gas stream (cm) pa is the aerosol particle density (g crnp3) and V is the particle velocity (cm s-'). (ii) Gravitational loss 4 = 3 (?,D/Pa@)+ (2) where dcg is the cut-off diameter (pm) of the aerosol droplet due to gravitational settling D is the chamber diameter (cm) g is the gravitational constant ( 10-3N m2 kgV2) and r is the time spent in the chamber. (iii) Turbulence loss dCt = ( 3 Q / W * ( 3 ) where d is the cut-off diameter (pm) of the aerosol droplet due to turbulence loss K is a constant whose value is deter- mined by the degree of turbulence of the system Lis the tube length and Q is the total gas volumetric flow rate.(iv) Centrifugal loss (4) where d is the cut-off diameter of the aerosol droplet due to the centrifugal loss M. is the inlet width (cm) n is the number of rotations and aerosol path taken (dimensionless) Vi is the inlet gas velocity and ps is the gas density (g ~ m - ~ ) . As the various interacting forces between the gas flow the aerosol stream and spray chamber are examined the variation of the aerosol flow velocity as the aerosol passes the chamber to the plasma should not be ignored. Since the Monte Carlo technique yields information on a 'model system' more useful and accurate results will be available if a more realistic model is assigned to the system.A jet model8 is thus proposed. A schematic diagram of the model for a round jet nozzle is shown in Fig. 1. The aerosol flow area is divided into moving and diffusion parts. Inside the former are the core zone where the aerosol flow velocity is equal and the mixed zone where the aerosol velocity attenuates and finally goes to zero at the boundary. The length of the core zone is determined by the shape of the nebulizer nozzle and is approximately equal to five times the diameter of the round jet nozzle (W,). In the diffusion area the aerosol velocity decreases as it moves down along the central part of the chamber and finally goes to zero at the boundary with increasing aerosol diffusion angle 2 is the axial distance and R is the radial distance (cm); Vo is the aerosol velocity (cm s - l ) at the jet nozzle i.e.at Z=O and V is the aerosol particle velocity at the 2 and R coordinate axes. Thus within the moving part V= V for the core zone and ZoR-W,(Zo-Z)/2 "[' -( 1.34W,Z Moving area * Diffusion area I I_ Fig. 1 Schematic diagram of the jet model for a round jet nozzle W width of jet nozzle; 2 distance of moving area; Z axial direction and R radial direction for the mixed zone. Within the diffusion part Assuming that the diameter of the jet nozzle is equal to the inner diameter of the carrier gas outlet and that the jet primary velocity of the liquid sample leaving the nozzle is equal to the outlet velocity of the carrier gas then Q Vo %+ (7) AJ3 where Q is the carrier gas flow rate (mls-') and A is the cross-sectional area of the gas outlet.According to this model the velocity of the particles at any position of the spray chamber can be calculated. By inserting the calculated velocity into eqn. ( l ) the cut-off diameter of the aerosol droplet due to impaction at any vectorially displaced position can be found. The cut-off diameter of the droplet for gravitational settling turbulence and centrifugal loss can also be found from their corresponding equations. It is important to note that the proposed jet model is based on a gross assumption; strictly speaking it is not appropriate for the enclosed jet within the spray chamber. However the results obtained indicate that the Monte Carlo simulation based on this jet model is valid. To perform the simulation it is necessary to know the distribution of the primary aerosol droplet size.According to Nukiyama and Tanasawa,' the primary aerosol droplet diam- eter is expressed as follows 585 (s)'.' [ q l [ 1:Q,]'-5 +597 - (8) D- - ,--I/ p (sp)O.' where D is the Sauter median diameter (pm) V is the velocity difference between gas and liquid flows to the nebulizer (m s- ') s is the surface tension of the liquid ( 1OP3N m-') p is the liquid density (gcmd3) ql is the liquid viscosity ( P ) and Q and Q are the volume flow rate of liquid and gas (crnp3 s)' respectively. The Sauter median diameter is defined as where D is the droplet diameter and n(D) is the number of droplets of diameter D and Do and D are the lower and upper limits of the distribution respectively. Generally Do is assumed to be zero. Thus after obtaining Ds from eqn.(8) the DM value can be calculated from eqn. (9). The primary aerosol droplet size distribution [O,D,] which is Gaussian can thus be obtained. In this work the parameters of three different concentric nebulizers with Scott-type spray chambers are used for simu- lation. The time interval for each simulation is 0.01 ms. For every Ar each particle is monitored for a vectorial displace- ment. Since the chamber is axially symmetrical the movement of the aerosol particle in each of the rectilinear coordinates is given by (11) where d and d are the axial and radial movement (cm) respectively; R and Rg2. are the random Gaussian number distributed about zero with a standard deviation of 1.0; V is the aerosol flow velocity (cm s-'); and D T is the temperature dependent diffusion coefficient (DT = D0(T/T0)3/2 where Do is the diffusion coefficient at To=273 K).The vectorial sum of these motions produces the new location of the particle. For every Af to determine whether the particle is lost by gravitational settling centrifugation or turbulence the normal random number R within the interval of 0 and D is compared dZ = Rgl (2D-&)+ + I/& dR = Rg2 (2D,At)* (12)JOURNAL OF ANALYTICAL ATOMIC SPECTROMETRY JUNE 1994 VOL. 9 703 with the corresponding cut-off diameter of the droplet due to the processes mentioned above. If R,3dc the particle is lost and drained to waste. If R,<dc the particle will then remain in the gas. From the cumulative movement it is possible to determine whether the particle collides with the chamber wall or with any other objects such as the paddle in the chamber.If impaction occurs the particle velocity can be calculated from eqns. (4) (5) and ( 6 ) depending on the position of the particle. Thence the cut-off diameter of the droplet dci due to the impaction loss can be obtained and compared with R,. Similarly if R 3 dci the particle is lost whereas if R < dci the particle will remain in the gas. When the movement Edz is greater than the chamber length the simulation is complete. Through the statistical cumulation the final particle loss number N1 and the particle number N2 that remain in the gas can be obtained. The mass transport efficiency E can be calculated from the folowing equation N N x 100% = A x 100% (13) N2 Nl+N2 En=- where N is the total particle number for simulation i.e.200 000 in this work. Hence the total analyte mass transport rate Kot can be derived from E provided that the Q1 is specified since Wot =~plQiC/100 where C is the analyte concentration (pg ml-I). r-1-1 Read in parameters Experimental The Monte Carlo simulation flow chart is depicted in Scheme 1. The program is written in FORTRAN 77 and is run on an M340 computer. It has a memory of 12 megabytes the hard disk has a memory of 4300 megabytes and the CPU compu- tation velocity is 2 500000 times s-'. To evaluate the validity of the program the simulation results are compared with those obtained experimentally by Browner et a1.' The dimensions of the nebulizer nozzles provided by Canals et aE." are used in the simulation work.They are listed in Table 1. The analyte solution used is 1000 pg ml-I Mn(NO,),. Results and Discussion Selection of Simulation Particle Number Generally the estimated value approaches the true value as the simulation particle number goes to infinity. However this would be too time consuming and is not realistic. Fig. 2 elucidates the relationship between the simulation particle number and the total mass transport rate. From the figure it can be seen that the total mass transport rate Wot becomes stable as the simulation particle number increases to 100000. In this work 200000 is selected as the particle number in order to guarantee higher precision. The total time required generated aerosol particles L v I Turbulence loss? increment + numberof loss particle t N,+N22N? Y I I N N I Increment of diffusion move (dZ dR)? N Increment number of particles to plasma N2 Parameters change? N End Scheme 1704 JOURNAL OF ANALYTICAL ATOMIC SPECTROMETRY JUNE 1994 VOL.9 Table 1 Dimensions of the nebulizer nozzles" Concentric nebulizer Parameter Inner diameter of the inner sample tube IJmm Outer diameter of the inner sample tube I,/mm Inner diameter of the outer gas tube IJmm Cross-section area of the gas outlet A,/mm2 Cross-section area of the liquid outlet AJmm' Recess of the inner tuber IJmm 1 0.318 0.510 0.580 0.060 0.079 -I- 0.07 - 2 0.508 0.700 0.730 0.034 0.203 '0.06 - 3 0.424 0.538 0.574 0.03 1 0.141 0.06 " 0.12 I I 0) 1. -. 5 0.11 : 0.10 g 0.09 r F CI " t E I I I I 1 10 20 30 40 50 Particle number x 10' Fig.2 Relationship between simulation particle number and total mass transport rate to obtain a simulation is in the range 8-25 min depending on the carrier gas flow rate used. Evaluation of the Validity of the Monte Carlo Program To evaluate the validity of our Monte Carlo program the following simulations were performed using the dimensions of the nebulizer nozzle listed in Table 1. Influence of Q on E and W, Simulation data of E and W, are obtained by increasing Q1 from 0.63 to 1.90 ml min-' while keeping the nebulizer gas flow rate Q constant at 0.65 1 min-'. The data are compared with the experimental data of Browner et aL9 and presented in Table 2. The results show that on the whole the simulation data agree quite well with those given in ref. 9. Influence of Q on WIot The value of Q was increased from 0.39 to 1.10ml min-' while keeping QI constant at 0.63 ml min- for three different concentric nebulizers (Nos.1-3). The simulation Wot data obtained are compared with the experimental results obtained by Browner et aL9 in Table 3. These results show that the simulated data correspond well with Browner's experimental results especially for lower values of Q,. Some discrepancy between them is observed for high values of Q but the trend that Fot increases with an increase of carrier gas flow rate is the same. Therefore it can be concluded that the model proposed is reliable and the program is valid. This program can be used to evaluate the performance of the nebulizer as well as to optimize experimental conditions.Aerosol Droplet Loss Mechanism Aerosol frow velocity distribution along the spray chamber In the jet model depicted in Fig. 1 it has been assumed that the aerosol flow velocity varies in the Z axial and R radial directions as the aerosol particles move along the spray chamber to the plasma torch. To illustrate this variation the following simulations are performed. The dimension of the No.1 nebulizer nozzle is taken into consideration. The values of Q and QI are kept at 0.65 and 0.63 ml min-' respectively. The results obtained by increasing the Z axial distance at R = 0 are presented in Table 4. The data indicate that as 2 increases from 0.1 to 14.5 cm the aerosol flow velocity decreases signifi- cantly from 180.6 to 3.6m s-'. The velocity decreases nearly 50% as Z increases from 0.5 to 1.0 cm and from 1.0 to 2.0 cm.To study the variation of the flow velocity along the radial distance simulations of the variation of the aerosol flow as a function of radial distance at various values of Z are performed. The results are presented in Table 5. The results show that generally the aerosol flow velocity V decreases with an increase of radial distance R at any axial distance 2; but varies significantly with R near the nebulizer nozzle. The flow velocity decreases almost 90% as the aerosol moves radially from 0.05 to 0.10 cm at Z=OS cm and from 0.10 to 0.20 cm at Z= 1 cm. However as 2 increases to 5 cm the relationship between V and R is no longer significant. The phenomena described above indicate that the proposed model might be helpful to give an understanding of the various loss processes taking place in the chamber.It can be presumed that it is the significant attenuation of the aerosol flow velocity that causes gravitational settling and turbulence and centrifu- gal and impaction loss processes to occur more easily. Table3 Effect of Q on W, Table 2 Effect of Q on E and W, Nebulizer QJI min-' Exp.* Simt 1 0.63 0.101 0.094 1.90 0.107 0.098 2 0.63 0.116 0.124 1.90 0.129 0.105 3 0.63 0.154 0.158 1.90 0.132 0.114 EXP Sim 0.97 0.90 0.39 0.3 1 1.1 1 1.18 0.4 1 0.33 1.52 1.51 0.43 0.36 Q 1 min 0.39 0.52 0.65 0.84 1.10 No 1. nebulizer ' Exp.* Sim.? 0.020 0.028 0.040 0.051 0.086 0.094 0.150 0.185 0.225 0.255 No 2. nebulizer Exp. Sim. 0.020 0.031 0.049 0.058 0.122 0.410 0.220 0.250 0.340 0.390 No 3.nebulizer Exp. Sim. 0.020 0.027 0.070 0.071 0.154 0.158 0.240 0.278 0.340 0.366 * Exp = experimental data from Browner et aL9 t Sim =simulation data from this work. * Exp = experimental data from Browner9 t Sim = simulation data from this work.JOURNAL OF ANALYTICAL ATOMIC SPECTROMETRY JUNE 1994 VOL. 9 705 Table 4 Variation of aerosol flow velocity as a function of 2 axial displacement Z/cm 0.1 0.2 0.5 1 .o 2.0 3.0 4.0 6.0 10.0 12.0 14.5 V/m s-' 180.6 180.6 104.7 52.4 26.2 17.5 13.1 8.7 5.3 4.4 3.6 Table5 Variation of aerosol flow rate velocity V as a function of radial distance R at various values of axial distance 2 Z/cm 0.5 1 5 10 Rlcm 0.05 0.10 0.13 0.14 0.05 0.10 0.20 0.25 0.27 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 vlm s-' 24.30 3.29 0.04 0 23.75 22.15 1.65 0.11 0 5.41 3.21 1.81 0.89 3.54 2.71 2.1 1 1.61 Effect of Qg on the Distribution of Loss Particles in the Chamber To determine which of the loss processes predominantly causes the separation of large droplets in the aerosol chamber the simulation of the number of lost aerosol particles Ni due to the various loss processes as a function of Q are performed.The dimensions of the No.1 nebulizer are used. The liquid flow rate Q1 is kept at 0.63ml min-l while the carrier gas flow rate is varied from 0.39 to 1.10 1 min-l. The total particle number Ntot for simulation is 200 000. The percentage of particles lost as a result of any loss process can be calculated from Ni/NtOt x 100. The data obtained are presented in Table 6. It is interesting to note that the percentage of particles lost and their distribution due to the various processes are influ- enced greatly by the magnitude of Q,.For example at lower Q the percentage loss due to turbulence processes is highest (87.65) but it decreases abruptly as Q increases to 0.84 1 min-l. On the other hand the percentage of particles lost due to impaction increases with an increase of Q,. At a Q of Table 6 Effect of Q on the distribution of loss particles Q,/1 min - 0.39 0.52 0.65 0.84 1.10 Q,/1 min-l 0.39 0.52 0.65 0.84 1.10 Centrifugal loss Turbulence loss No. YO 15 520 7.76 1169 0.58 0 0 0 0 0 0 Gravitational loss No. YO 175 290 87.65 164 861 82.43 109 000 54.50 21 693 10.85 785 0.39 Impaction loss No. Yo 2 770 1.39 16 061 8.03 49 062 24.53 75 107 37.55 21 720 10.86 No. % 6 000 3.00 17 909 8.47 41 938 20.07 103 200 49.84 177 495 86.32 1.10 1 min-' the percentage loss reaches 86.32. Gravitational settling also plays a role in the separation of large droplets. At Q values of 0.65 and 0.841 rnin-l the percentage losses due to gravitational settling are 24.53 and 37.55 respectively.The role of centrifugal loss is insignificant. At low Q the percentage loss due to centrifugation is only 7.7 and it goes to zero with an increase of Q,. Evidently the predominanting loss process in the separation of the large droplets in the aerosol chamber depends strongly on the magnitude of Q,. When a lower Q is applied i.e. from 0.39 to 0.52 1 min-l it appears that the turbulence-induced loss is the predominant mechanism whereas when a higher Q is applied e.g.1.1 1 min - ' the impaction-induced loss might predominate. Conclusions A jet model is proposed for calculating the variation of the aerosol flow velocity in a spray chamber. The Monte Carlo program is evaluated by comparing the simulated and exper- imental data of E and KO,. Good agreement between them is obtained. This suggests that the proposed model is reliable. The program can be used to select and evaluate the perform- ance of nebulizers and to optimize the operating conditions for ICP-AES work. The aerosol flow velocity distribution is a function of the axial and radial distances of the chamber. The results indicate that the proposed velocity model might be useful in providing an understanding of the various loss processes that occur in the chamber.The effect of the carrier gas flow rate on the distribution of loss particles due to various loss processes along the chamber has been studied. The results indicate that the predominant loss mechanism depends strongly on the magnitude of Q,. However the results obtained are only preliminary. In order to obtain the relation between signal intensity and aerosol particle behaviour the vaporization atomization and ionization processes in the ICP have been studied and simulated. The results will be reported in a future paper. This work was supported by the National Natural Science Foundation of China. References 1 2 3 4 5 6 7 8 9 10 Nukiyama S. and Tanasawa Y. Trans. SOC. Mech. Eng. (Jpn.) 1939 4 68. Browner R. F. Boorn A. W. and Smith D. D. Anal. Chem. 1982,54 1411. Skogerboe R. K. and Olson K. W. Appl. Spectrosc. 1978,32,181. Gustavason A. Spectrochim. Acta Part B 1984,39 85. Sharp B. L. J. Anal. At. Spectrom. 1988 3 613. Sharp B. L. J. Anal. At. Spectrom. 1988 3 939. Zheng J. Zhang Z. and Qian H. Proceedings of the Fourth International Beijing Conference and Exhibition on Instrumental Analysis C. Spectroscopy 1991 p. 13. Seichi H. and Yotab O. Jet Engineering Science Publishing House China 1977 66. Browner R. F. and Canals A. and Hernandis V. Spectrochim. Acta Part B 1992 47 659. Canals A. Hernandis V. and Browner R. F. Spectrochim. Acta Part B. 1990 45 591. Paper 3103951 J Received July 7 1993 Accepted January 12 1994