10524
Ind. Eng. Chem. Res. 2010, 49, 10524–10534
Evaluation of Photon Absorption in an Aqueous TiO2 Slurry Reactor Using Monte Carlo Simulations and Macroscopic Balance Jesus Moreira, Benito Serrano,† Aaron Ortiz, and Hugo de Lasa* Faculty of Engineering Science, Chemical Reactor Engineering Centre, The UniVersity of Western Ontario, London Ontario, Canada. N6A 5B9
The radiation field in an annular photocatalytic reactor is simulated using a Monte Carlo method (MC) for two TiO2 suspensions in water. Simulations are performed by using both the spectral distribution and the wavelength-averaged scattering and absorption coefficients. The Henyey-Greenstein phase function is adopted to represent forward, isotropic, and backward scattering modes. It is assumed that the UV lamp reflects the backscattered photons by the slurred medium. Photoabsorption rates using MC simulations and spectral distribution of the optical coefficients agree closely with experimental observations from a macroscopic balance. It is found that the scattering mode of the probability density function is not a critical factor for a consistent representation of the radiation field. MC simulation for the optimal catalyst concentration reveals that the maximum LVREA is reached at a concentration of 0.14 g L-1 for TiO2 Degussa P25. From this concentration, the apparent optical thickness is determined to be 2.8476 which is in agreement with the optimal one previously reported. This concentration is comparable to that determined experimentally for phenol photocatalytic degradation. 1. Introduction There is a wide variety of technical fields were TiO2 heterogeneous photocatalysis is of importance. Among its applications, self-cleaning surfaces, water purification, air purification, self-sterilizing surfaces, antifogging surfaces, anticorrosion applications, environmentally friendly surface treatment, photocatalytic lithography, etc. can be mentioned.1,2 For water decontamination, heterogeneous photocatalysis is considered a tertiary water treatment process that belongs to the family of advanced oxidation process (AOP) technologies. This new technology is a promising alternative method for the removal of organic and inorganic pollutants in water.3 Photocatalytic reactions are the result of the interaction of light of appropriate energy with a solid semiconductor.4 When the incident photon has energy equal or greater than the semiconductor band gap, radiation is absorbed and electrons are moved from the valence band to the conduction band giving rise to the formation of electron-hole pairs.5 The resulting electron hole pairs promote redox processes that may completely mineralize and transform organic compounds into water, carbon dioxide, and some mineral acids. Among the TiO2 semiconductor photocatalysts tested, Degussa P25 (DP25) has proven to be one of the most active catalyst with very interesting attributes such as high stability, good performance, and low cost.6 Slurry photoreactor systems have shown the largest photocatalytic activity when compared to photocatalytic reactors with immobilized catalyst.3 Advantages of slurry photoreactors include fairly uniform catalyst distribution, high photocatalytic surface area to reactor volume ration, limited mass transfer, well-mixed particle suspension, and low pressure drop throughout the reactor. The main disadvantage for slurry photoreactors is the requirement of postprocess filtration for catalyst recovery/recycling.3,6 Efficiency determinations in photocatalytic reactors allow for a comparison between experimental results obtained from * To whom correspondence should be addressed. E-mail: hdealsa@ eng.uwo.ca. Phone: 519 661 2111ext 82144. Fax: (519) 661 3498. † Present address: Unidad Acade´mica de Ciencias Quı´micas, Programa de Ingenierı´a Quı´mica, Universidad Auto´noma de Zacatecas, Zacatecas, Zac., Me´xico.
different laboratories and under different experimental conditions.7 Several efficiency definitions are found in the technical literature. The most frequent parameter is the quantum yield,8 which relates the radicals produced on the catalyst surface during the primary reaction processes per absorbed photon. Serrano and de Lasa9 also proposed a photochemical thermodynamic efficiency factor (PTEF). This parameter relates the energy needed to produce OH• radicals over the irradiated energy absorbed by the photocatalyst. In either case, the efficiency determination involves the same key variable, the rate of photons absorbed by the photocatalyst. Photoreaction can take place only if the light involved in the process has enough energy to promote TiO2 excitation (3.2 eV for anatase). Therefore, an accurate estimation of light intensity distribution is critical in the design and rating of photoreactors.10 This task involves the estimation of the local volumetric rate of energy absorption (LVREA), which is derived from the radiative transfer equation solution (RTE).11 The RTE is an integro-differential equation that describes the light intensity distribution in a photoreactor. For homogeneous photoreaction systems, an analytical solution for the RTE is feasible.12 For a heterogeneous medium however, and given that solid particles cause scattering and light absorption, an analytical solution is only possible under simplified assumptions.13-15 Furthermore, due to the heterogeneous nature of TiO2 particles, scattering occurs according to mechanisms that are quite different from those in multiphase gas-liquid systems. In photocatalytic systems, the light intensity distribution in a photoreactor is a function of: (a) lamp type, (b) reactor geometry, (c) type of catalyst, (d) catalyst concentration, (e) particle agglomeration size, (f) the nature of the reactor walls (highly reflecting reactor walls, specular and diffuse reflecting walls, nonreflecting walls), (g) flow rate, (h) pH, (i) recycle flow rate, and (j) radiation wavelength.7,16 Therefore, a numerical solution of the RTE seems to be a viable alternative. The radiation intensity inside a photoreactor can also be determined experimentally. A simple method is chemical actinometry, which gives high values for the LVREA because it assumes that the reacting medium absorbs all photons reaching the inner reactor wall. Salaices et al.7,17 proposed an alternative
10.1021/ie100374f 2010 American Chemical Society Published on Web 05/12/2010
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
experimental method for the evaluation of the rate of photon absorption inside an annular photoreactor. This method uses radiometric measurements along with a macroscopic radiation balance to determine the radiation absorbed by the solid catalyst, as well as the forward and backward scattering radiation and the extinction coefficients. The advantage of this method is that the evaluation of the optimal catalyst concentration can be determined by a single spectrometric measurement of the forward radiation outside the reactor (0.2% transmission of the total input radiation). This experimental method, however, needs to be validated with a rigorous numerical calculation of the LVREA. The LVREA distribution can be determined in any computational domain by solving the RTE. Several issues have to be addressed before this can be accomplished. First, the optical properties of the reaction medium (i.e., absorption and scattering coefficients) have to be known.10 Second, the boundary conditions (light being received by the radiation source) have to be precisely established. The most common numerical methods for finding the solution to the RTE are the discrete ordinate method (DO), the finite volume method (FV), and the Monte Carlo method (MC). Stochastic simulation methods, such as the MC method, are preferred over deterministic methods for finding the LVREA for complicated geometries.18 Yokota et al.12 states that a statistical approach to assessing the absorbing and scattering phenomena in heterogeneous systems is the most effective tool. In MC simulations, individual photons (or bundles of photons) are traced from their creation until the photons are either absorbed or scattered from the system. Pareek et al.10 presented an MC approach to estimate LVREA in an annular photoreactor. These authors divided the reaction space into small cubic cells and assumed that the lamp surface is a uniform line source emission. Changrani et al.18 employed an MC methodology to simulate a 3D polychromatic UV radiation field in an annular photocatalytic reactor. Changrani et al. used alumina foams as the catalyst support structure. Yokota et al.12 developed an MC simulation model having three parameters: an attenuation coefficient, a probability of ray absorption, and a scattering parameter for isotropic and anisotropic scattering modes. Yokota et al. used this model to predict the light intensity in a photoreactor with heterogeneous medium; they found that simulation results were close to the observed experimental data. In the present study, MC simulations are performed to determine the LVREA for the Photo-CREC Water-II reactor with two different types of TiO2 photocatalysts (Aldrich and Degussa P25). An experimental macroscopic balance is used to determine the light distribution inside the reactor. In the simulation, four parameters were used: (i) an extinction coefficient, (ii) the probability of photon absorption, (iii) the probability of backward scattering absorption by the internal Pyrex glass tube, and (iv) the Henyey-Greenstein phase function describing the forward, isotropic, and backward scattering. It is also assumed that the inner UV lamp reflects photons backscattered by the TiO2 suspension. Comparisons of light absorption and forward scattering radiation rates from experimental macroscopic balances and MC simulation show a satisfactory agreement. These findings demonstrate the effectiveness of the MC simulation for the assessment of photon absorption in the photocatalytic reactor and, also, its usefulness for determining the optimal catalyst concentration from the macroscopic radiation balance.
10525
2. Radiation Transport Equation The application of the RTE to photocatalytic processes is done by making a radiation balance across a thin slab. The resulting equation may be expressed as16,19,20 dIλ(s, Ω) ) -κλIλ(s, Ω) - σλIλ(s, Ω) + j eλ(s) + ds 1 σ p(Ω′ f Ω)Iλ(s, Ω′) dΩ′ 4π λ 4π
∫
(1)
where Iλ is the spectral specific intensity of radiation having a wavelength λ (einstein m-2 s-1 sr-1), κλ is the absorption coefficient (m-1), σλ is the scattering coefficient (m-1), and p(Ω′ f Ω) is the phase function for the in-scattering of photons. The first term on the right-hand side of eq 1 represents the loss of photons due to absorption, the second term considers the loss of radiation due to out-scattering, the third term accounts for emission of light due to temperature effects, and the fourth term shows the gain in radiation due to in-scattering. Integration of this partial-integro-differential equation requires at least one boundary condition at the point of radiation entrance to the reactor volume which can be provided by an appropriate lamp emission model. Pareek et al.10 presented a good review of line source, surface source, and volume source models. They observed that, in general, the volume source model gives the best results for different situations, but line source models may give reasonably accurate results for annular photoreactors. Due to photocatalysis being carried out at low temperatures, the emission term can be neglected. Thus, the sum of κλ + σλ is called the extinction coefficient βλ. Equation 1 could therefore be rearranged to give dIλ(s, Ω) 1 σ ) -βλIλ(s, Ω) + ds 4π λ
∫
4π
p(Ω′ f Ω)Iλ(s, Ω′) dΩ′
(2) If the local incident radiation at any point from all the directions is given by Gλ(x, y, z) )
∫
Ω)4π
Ω)0
Iλ(s, Ω) dΩ
(3)
Then, the LVREA at any point can be represented as El,λ(x, y, z) ) κλGλ(x, y, z)
(4)
The absorption threshold for TiO2 depends on the energy band gap. For DP25, the band gap is 3.2 eV, and those photons with a wavelength less than or equal to the band gap energy (λEbg ) 390 nm) promote excitation of electrons in the semiconductor particles. Therefore, the total LVREA is given by El )
∑
λeλbg
El,λ(x, y, z) )
∑
κλGλ(x, y, z)
(5)
λeλbg
Since the values for the κλ and σλ coefficients for titania powders depend on the wavelength of the light source, the RTE has to be solved for each individual wavelength. In most cases, this complicates the solution of the RTE, not only because it requires extensive use of computer memory, but also because the determination of absorption and scattering coefficients is not a trivial exercise. Romero et al.21 presented the wavelengthaverage values for the absorption and scattering coefficients of DP25 and Aldrich TiO2 in the 0.05-0.5 g L-1 concentration range for a Hanovia LL-189a-10/1200 lamp with an emission range of 220-1370 nm. For this study, a Black Blue Light Lamp
10526
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
with an emission range from 315 to 420 nm is used (see Figure 3), the averaged values for the absorption and extinction coefficients were determined from the data reported by Romero et al.22 The method used to obtain these coefficients is described by Toepfer et al.23 These values are respectively σDP25 ) 5.3707Wcat κDP25 ) 0.6748Wcat
(6a)
σAldrich ) 3.2977Wcat κAldrich ) 0.5985Wcat
(6b)
with σi and κi being the wavelength-average scattering and absorption coefficients in inverse meters for the respective titania powder, and Wcat as the catalyst loading in grams per cubic meters. The specific absorption coefficients for Aldrich and Degussa P25 as a function of wavelength (315-435 nm) were reported by Romero et al.22 These discrete values, along with wavelength-average scattering and absorption coefficients are used in this study, for the determination of LVREA by the MC method. Numerical simulations are compared with experimental observations. In photocatalysis, there is multiple scattering involved due to the topography of TiO2 particles. The parameter describing the scattering mode in eq 2 is the phase function p(Ω′ f Ω). This parameter gives the probability that a photon will be scattered from the direction Ω′ to the direction Ω. Therefore, the selection of the phase function is an important step in any calculation where multiple scattering is involved. Computing the new directions of the scattered photons is perhaps the most challenging task in MC simulations, requiring a large amount of computer time.24 Thus, complicated phase functions require a large computation time leading to inefficient simulations. In an established scattering problem however, the phase function is given, not chosen. Therefore, it is customary to use a phase function that preserves the main characteristics of the real phase function while at the same time rendering manageable computation of the scattering angles.25 For this purpose, the Henyey-Greenstein (HG) phase function seems to be an appropriate choice. This is a one parameter function that is able to reproduce a wide range of scattering probability density functions.26
Figure 1. Schematic representation of the Photo-CREC Water-II reactor.
Figure 2. Dimensions for the Photo-CREC Water-II reactor.
3. Experimental Section A StellarNet EPP2000C-25 LT16 spectrometer was used to perform spectroscopic measurements in the Photo-CREC WaterII reactor. This reactor was previously described by Ortiz-Gomez et al.27 A schematic representation of the Photo-CREC WaterII reactor is reported in Figure 1. The reactor is constituted by the following components: (1) black light lamp or (BL lamp), (2) Pyrex glass inner tube with diameter of 3.58 cm, (3) replaceable Pyrex inner tube with diameter of 5.6 cm, (4) silica windows, (5) black polyethylene outer tube, (6) stirred tank, (7) centrifugal pump, (8) air injector, and (9) sampling port. Seven 1.1 cm in diameter circular windows made of fused silica are equally spaced along the reactor outer cylinder wall to allow radiation transmission measurements. The external cylinder, represented by the number 5 in Figure 1, was made of nonreflecting black polyethylene in order to eliminate the radiation reaching the inner surface wall that could be reflected toward the sensor cell. Also, the Photo-CREC Water-II reactor has nonreflecting walls, and this means that every time a photon of light reaches the outer reactor wall, the photons are absorbed and its trajectory is terminated. The dimensions of the photo-
Table 1. Dimensions for the Photocatalytic Reactor and Lamp Characteristics component annular reactor
black blue light lamp (UVP-XX-15BLB)
parameter
value
internal radius external radius height internal Pyrex glass thickness input power
1.76 4.44 44.5 0.23
cm cm cm cm
output power length radius emission range emission rate
4W 41.3 cm 1.33 cm 300-420 nm 1.1910 × 10-5 einsteins s-1
15 W
reactor presented in Figure 1 are now reported in Figure 2. These values are important for the MC simulation. Finally, the characteristics for the lamp used in this study are presented in Table 1 along with a summary for the dimensions of the reactor. Titanium dioxide suspensions for radiation transmission measurements in the Photo-CREC Water-II were prepared with distilled water. Before any measurement, the reactor was
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
10527
Table 2. Lamp Emission and Pyrex Glass Transparency Rates for the BL Lamp variable Po Pi Pa-wall
einsteins s-1 -5
1.1910 × 10 1.1210 × 10-5 6.9995 × 10-7
Watts 4.034 3.778 0.246
(c) Pbs is approximated from the difference between Pi and the rate of photons transmitted when the catalyst concentration approaches zero (assumptions and explanations for this approximation are reported in Salaices et al.7 and references therein): Figure 3. Typical radiative flux spectra for a black light lamp. The full line represents the emission spectra for the lamp; the dashed line represents the radiation spectra variation through the 0.23 cm-thickness × 3.05 cmdiameter Pyrex inner tube.
thoroughly washed to ensure that no foreign particles were present during any measurement. The radiation transmission was measured first for an empty reactor with and without internal tube in order to measure the transparency of the internal Pyrex tube. The following step is to fill out the reactor with 6 L of distilled water and then increase the TiO2 concentration from 0 to 0.2 g L-1. For every catalyst concentration, the radiation transmission in the seven windows (see Figure 1) was done. When building the radiation transmission profiles, measurements for windows 2-6 were averaged. Windows 1 and 7 were not considered because the radiation transmission was low. 4. Radiation Transmission Modeling The following sections present the macroscopic balance and the Monte Carlo approach used to determine the radiation absorbed by the solid photocatalyst inside the Photo-CREC Water-II. Macroscopic Radiant Energy Balance. A macroscopic radiant energy balance applied to a control volume in the photoreactor with boundaries containing the slurred TiO2 catalyst (as proposed by Salaices et al.7) was performed to estimate the rate of photon absorption: Pa ) Pi - Pbs - Pt
(7)
where Pa is the rate of absorbed photons, Pi is the rate of photons reaching the reactor inner surface, Pbs is the rate of backscattered photons exiting the system, and Pt is the rate of transmitted photons in einsteins per second. The various terms in eq 7 are estimated as follows: (a) Pi is estimated from the rate of photons emitted by the black light lamp (Po) minus the rate of photons absorbed by the inner Pyrex glass (Pa-wall): Pi ) Po - Pa-wall
(8)
The rate of photons emitted by the lamp can be determined from radiometric measurements and from the lamp emission spectrum (see Figure 3). Po )
∫ λ∫ ∫ λ2
λ1
L
0
2π
0
q(θ, z, λ)r dθ dz dλ
(9)
(b) Pa-wall is estimated from transmission measurements through the inner Pyrex tube. Figure 3 reports the change of spectrum emission as radiation evolves through the wall of the inner Pyrex glass. The location of the λEbg is also presented.
Pbs ) Pi - P| Cf0+
(10)
(d) Pt is the addition of the transmitted nonscattered radiation (Pns) and the forward-scattering radiation (Pfs). Pt ) Pns + Pfs
(11)
One should mention that Pfs and Pns can be estimated by using polished-aluminum tube collimators, which account for the combined transmitted nonscattered radiation and forwardscattering radiation, and black collimators, which account only for the transmitted nonscattered radiation, respectively. Monte Carlo Simulations. The radiation scattering mechanism in a photocatalytic reactor starts when a photon is emitted from the lamp, travels a distance l, and then is either absorbed or scattered within the reacting medium. The generated photons interact with the reacting medium according to probabilistic interactions determined by the scattering and absorption cross sectional areas (σλ and κλ) of the reacting medium and also the phase function. The MC simulation begins with a given total energy input which is a function of the lamp used in the photocatalytic process. The spectrum intensity of the lamp used in this study is reported in Figure 3. Table 2 shows the lamp emission rate and the emission rate reaching the reactor inner Pyrex glass. It can be seen that 6% of the emitted radiation by the lamp is absorbed by the inner Pyrex glass. In the MC simulation, the RTE is solved by tracing trajectories and fates of all the photons emitted by the lamp. These trajectories and fates are found by using random numbers.12,16 Different MC simulations for radiation modeling could be considered depending on the underlying hypotheses. On the basis of the considerations given above the following model assumptions are adopted in this study: (i) Emission of photons from the lamp surface is assumed to be uniform and directionally independent. (ii) Emission of photons from the lamp surface is considered to be a stochastic process. (iii) Photons emitted by the lamp have a defined probability to be absorbed by the inner Pyrex glass tube before entering the reacting medium. This probability is defined by the transparency of the Pyrex tube (see Figure 3 and Table 2). The photons reflected to the lamp by the slurry medium have the same chances of being absorbed by the inner Pyrex tube. (iv) Photons scattered by the reacting medium are determined by the HG phase function. Forward, isotropic, and backward scattering are considered for the simulation, with photon reflection being assumed elastic. (v) Photons reaching the outer polyethylene reactor wall are considered absorbed by the reactor wall. As a result, these photons are counted as transmitted (i.e., its trajectory is terminated at the wall).
10528
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
2) U(1.6 cm, L + 1.6 cm) while R2 (R2 < 1) fixes its circumferential position U(0, 2π). (2) Once the photon emission coordinates are set, the direction of the photon flight is established in spherical coordinates by two angular coordinates; the zenith (θ) and azimuth (φ) angles, by using the HG phase function: pHG(θ) )
1 - g2 1 2 4π (1 + g - 2g cos(θ))3/2
(12)
where θ is the scattering angle and g is the asymmetry factor of the scattered radiation distribution. The g parameter of eq 12 can be set in the range of -1 < g < +1 considering in this manner phase functions ranging from a completely backward to a completely forward form. When g ) 0, eq 12 represents an isotropic phase function. The random event for which the zenith angle falls (with a probability density function given by eq 12) for the [θ, θ + dθ] interval24 is calculated by using a random number R uniformly distributed in the range [0, 1] such that
∫
θ
Figure 4. 3D view of the annular region used for MC simulation.
Assuming that a photon emitted in a random direction from the lamp surface enters the annular region of the reactor (see Figure 4), this photon travels a distance l (l g 0) with a finite probability of arriving to a point “A” through the heterogeneous medium without being scattered or absorbed. Once the photon reaches point A, there are two possibilities for the light ray:17 (a) the photon is absorbed and thus, its course is arrested; (b) it is scattered according to the HG phase function and its flight continues until this photon is either absorbed by a catalyst particle or reaches the reactor wall and its trajectory is terminated. On the other hand, if the photon is back-reflected toward the lamp, after traveling a distance l, the photon can be re-emitted at the same axial position but with different equatorial and angular angles. Mathematical Steps Involved. For this simulation, experimental results were used to model the spectral distribution of the BL lamp. The number of photons associated with each wavelength was experimentally determined for the range 260 nm e λ e 440 nm using spectrometric measurements. The rate of photon emission was also established at 7.2 × 1018 photons per second from spectrometric measurements. This experimental data was used as a starting point for the simulations. Once the number of photons was determined for every wavelength, their fate was traced using the MC simulations. In this way, the number of events considered in the MC calculations becomes dependent on the number of photons at every wavelength. Simulations were performed by using the averaged wavelength for absorption and scattering coefficients (eq 6) and also the spectral distribution of such coefficients as reported by the literature.22 The annular region presented in Figure 4 was divided into small cubic cells. Every time a photon was absorbed in a cubic cell, its value was saved in the corresponding cubic cell so that the LVREA could be calculated. The whole annular region was considered for the MC simulation, using a coordinate system as described in Figure 4. The trajectories of the photons for every wavelength were traced using the following steps in the context of the MC simulation: (1) The photon emission location on the lamp surface is determined first by using two random numbers: R1 (R1 < 1) sets the location along the z coordinate (refer to Figure
o
pHG(θ′) dθ′ ) R
(13)
In order to calculate the zenith angle, eq 13 has to be solved to obtain a solution that expresses θ as a function of R, with this equation being solved numerically. However, this approach increases considerably the computation time in MC simulations. As an alternative, a probability density function can be conveniently found by slightly modifying eq 12 as24 p˜HG(cos(θ)) )
1 - g2 1 2 (1 + g2 - 2g cos(θ))3/2
(14)
Therefore, a distribution for p˜HG(cos(θ)) can be represented as
∫
cos(θ)
-1
p˜HG(cos(θ′)) d(cos(θ′)) ) R
(15)
With eq 15 having an exact analytical solution expressed by
{[
(
1 1 - g2 1 + g2 cos(θ) ) 2g 1 - g + 2gR 2R - 1,
)] 2
, if g * 0 if g ) 0 (16)
Two random numbers (R3 and R4) uniformly distributed in the range [0, 1] are generated to calculate the zenith and azimuth angles. When cos(θ) is calculated with eq 16 by using R3, the zenith (latitude) angle for the photon flight direction is obtained by computing the θ ) a cos((cos(θ)). For the azimuth angle, the same probability of reflection is assigned; therefore, the scattering angles are calculated as follows: θ ) a cos(cos(θ)) φ ) 2πR4
(17)
In the simulations, different values for g were studied to include the different scattering modes. For isotropic scattering, g ) 0 and eq 17 becomes θ ) a cos(2R3 - 1) φ ) 2πR4
(18)
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
(3) Photons have a 6% probability of being absorbed by the inner Pyrex glass tube. This probability is determined by the transparency of the tube material (refer to Figure 3 and Table 2). If the photon is absorbed by the Pyrex glass, its trajectory is arrested and a new photon is generated (step 1). If at any time, the photon of light is backscattered by the suspension toward the lamp, the photons again have a 6% probability of being absorbed by the Pyrex glass. If the photon is not absorbed by the inner Pyrex glass tube, it is considered lamp reflected at the same axial position but with a different angle. (4) Once the photon emitted in the first steps enters the reaction zone, it travels a distance l without an interaction occurring along this path. Then, the next step is the evaluation of the photon flight length l. The probability of this event is given by10 P(l) ) e-βλl
(19)
where βλ is the extinction coefficient of the medium. Therefore, the flight length l can be calculated by a random number R5 (R5 < 1) as follows: l)-
1 ln(R5) βλ
(20)
Yokota et al.12 present a slightly different definition for the free path length, l)-
1 ln(1 - R5) βλ
(21)
This definition, however, renders the same results for the simulation because (1 - R5) is also a random number between 0 and 1. If after traveling a distance l, the photon is located at point A, its Cartesian coordinates are simply determined by
xnew ) xold + exl ynew ) yold + eyl znew ) zold + ezl
10529
(22)
where old refers to the previous location of the photons inside the reactor and new refers to the new location once the photon travels a distance l. The utilization of Cartesian coordinates for establishing photon location presents computational advantages because a photon’s direction is uniquely specified by the direction cosines (ex, ey, ez) with respect to the coordinate axes.18 (5) If the position of the photon after traveling a distance l is inside the annular region, two possibilities can happen: either the photon is absorbed by the medium or it is scattered to a new location. This step involves the probability calculation of the photon being absorbed. For photons crossing in the reactor annulus, their fate is determined by an absorption criterion, which is the probability that the photon is scattered. This absorption criterion is given by18 P(a) )
κλ κλ ) βλ κλ + σλ
(23)
Thus, another random number is generated (R6 < 1). If P(a) > R6, then the photon is absorbed and stored in the corresponding volume cell. At this point, the photon trajectory is terminated and the sequence of calculations is reinitiated for a new photon emitted by the lap surface (step 1). Otherwise, the photon is scattered and a new direction for the photon is established (step 2). However, if the photon of light is outside the outer polyethylene tube the photon is allowed to escape as a forward scattered photon and is assumed to be absorbed by the nonreflecting wall. In the MC simulations a large number of events must to be considered until the physical properties under investigation have
Figure 5. Experimental results for the transmitted radiation and the LVREA and comparison with “MC simulation 1” when the wavelength-averaged scattering and absorption coefficients are used and isotropic scattering is assumed: (O) experimental results, (s) MC simulation; (i) LVREA for DP25, (ii) Pt for DP25, (iii) LVREA for Aldrich, and (iv) Pt for Aldrich.
10530
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
Figure 6. Experimental results for the transmitted radiation and the LVREA and comparison with “MC simulation 2” when the wavelength-averaged scattering and absorption coefficients are used: (O) experimental results, (s) MC simulation; (i) LVREA for DP25, (ii) Pt for DP25, (iii) LVREA for Aldrich, and (iv) Pt for Aldrich.
Figure 7. Experimental results for the transmitted radiation and the LVREA and the comparison with “MC simulation 3” when the spectral distribution for the scattering and absorption coefficients is considered: (O) experimental results, (s) MC simulation; (i) LVREA for DP25, (ii) Pt for DP25, (iii) LVREA for Aldrich, and (iv) Pt for Aldrich.
small statistical fluctuations. Ideally, the number of events traced should be equal to the number of photons emitted by the lamp per unit time.10 However, this is a demanding computational process because the BL lamp emits of 7.2 × 1018 photons per second (1.19 × 10-5 einstein s-1). In practice, fewer events can be used for accurate calculations. Yokota et al.,12 for example, considered 105 events in an MC simulation for the prediction of light
intensity decay in an heterogeneous medium. Pareek et al.10 grouped the photons in packets so that the total number of events could be simulated in a personal computer. Using a similar approach, for the simulations in this study, 7.2 × 1018 photons were accounted forming 107 packets of photons, and as a result, events could be calculated using an Intel Core Duo PC (2 GHz). MC simulations of photon interactions employ random numbers to choose points of emission, optical depths, scattering
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010 Table 3. Least Square Error Calculation for MC Simulations When Isotropic Scattering Is Considered DP25 case
LVREA error -7
Aldrich Pt error
LVREA error
Pt error
-7
simulation 1 1.1698 × 10 2.2821 × 10 7.9346 × 10 1.5893 × 106 simulation 2 1.9780 × 10-5 1.9332 × 106 1.6892 × 10-5 1.4674 × 105 simulation 3 2.0850 × 10-8 5.2536 × 104 6.1326 × 10-8 4.1121 × 105 5
angles, absorption, and scattering probabilities. Since all these random numbers are generated by a set of algorithms in a computer, no output is truly random. Therefore, in order to produce sequences of numbers that pass a suitable randomness test, an algorithm has to be developed. The RAND function in MATLAB provides an excellent and easy way to generate pseudorandom numbers for MC simulations.10 This function (with a period of (219937 - 1)/2), exceeds any of the computational number of simulation steps. Therefore, the RAND number was considered suitable for MC simulations of the LVREA. Thus, for the simulations presented in this study, MATLAB programs were developed to estimate the absorbed photons in the Photo-CREC Water-II reactor with the RAND function being used for the generation of random numbers in all cases. 5. Results and Discussion Various parameters involved in the macroscopic radiation balance were determined experimentally based on the method proposed by Salaices et al.7 In this regard, the forward-scattered radiation at different titanium dioxide concentrations was measured using small tubes with a dark inner surface located at various axial reactor locations. Photon absorption and the LVREA were then evaluated for both Aldrich and DP25 catalysts for a 0-0.2 g L-1 concentration range. One proposed MC method simulation (simulation 1) to predict the LVREA assumed that those photons that are backscattered by the medium and impinge on the lamp are reflected at the same axial position, (i.e., there is no absorption of photons by the fluorescent BL
10531
lamp). On the other hand another MC simulation assumed that photons reflected by the medium impinging on the lamp were actually absorbed by the lamp (simulation 2). For the abovedescribed simulations, the wavelength-averaged scattering and absorption coefficients were used (eq 6) and isotropic scattering was assumed. Equation 18 was employed to calculate the reflecting angles for the photons inside the reactor (i.e., g ) 0). Radiation absorption experimental results and MC simulations (simulation 1) are reported in Figure 5. MC results for simulation 2 are reported in Figure 6. The LVREA predictions presented in Figure 5 confirm that the use of wavelength-averaged parameters (eq 6a) renders good prediction for the LVREA. However, the prediction for the LVREA is not very accurate, in Figure 5, one can notice an overestimation for the LVREA. Given that “simulation 1” provides a more realistic scenario, it was chosen for the evaluation of the LVREA in the PhotoCREC Water-II reactor. However and to assess the influence of the scattering and absorption coefficients, a “simulation 3” was performed. Figure 7 reports these results when the spectral distribution for the scattering and absorption coefficient is considered with isotropic scattering (g ) 0) assumed. Comparing Figures 5 and 7, one can conclude that by using the spectral distribution for the MC simulations, a more accurate prediction for the LVREA and the transmitted radiation inside the reactor is obtained. Table 3 summarizes the errors from experimental data and MC simulations 1, 2, and, 3 adopting the isotropic scattering assumption. Knowing that simulation 3, involving spectral distribution of scattering and absorption coefficients, renders the smallest errors, simulation 3 is proposed for further calculations to establish the significance of asymmetry factors. With this end “simulation 4”, involving a variation of the asymmetry factor was implemented with this factor changing from a narrow forward peak (g ) 1) to a narrow backward peak (g ) -1). Results for different values of g are shown in Figure
Figure 8. Asymmetry factor influence in MC simulation 4: (O) experimental results, (- - -) g ) 1, (- - -) g ) -1, (green s) g ) 0.8 and -0.8, (blue s) g ) 0.5 and -0.5, and (red s) g ) 0.1 and -0.1; (i) LVREA for DP25, (ii) Pt for DP25, (iii) LVREA for Aldrich, and (iv) Pt for Aldrich.
10532
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
Table 4. Least Square Error Calculation for MC Simulations for Different Asymmetry Factor Values Error in MC Simulation 4 from Experimental Values DP25 g 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 isotropic -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
LVREA error 1.8290 5.9836 3.5676 1.9869 1.1765 7.1290 4.7211 3.4762 2.9021 2.4506 1.8044 1.9001 2.0462 2.0553 3.1621 5.1632 8.9136 1.5784 2.8901 4.8954 1.0594
× × × × × × × × × × × × × × × × × × × × ×
10-6 10-7 10-7 10-7 10-7 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-7 10-7 10-7 10-6
Aldrich Pt error
7.8864 4.4372 2.5181 1.4621 9.6012 7.3418 6.4109 6.4542 6.6325 6.7952 5.2536 7.4769 7.8006 8.1381 8.7256 9.5548 1.1239 1.5392 2.3516 3.9029 7.5337
× × × × × × × × × × × × × × × × × × × × ×
105 105 105 105 104 104 104 104 104 104 104 104 104 104 104 104 105 105 105 105 105
LVREA error 3.5191 6.1424 3.2302 1.7964 1.0643 7.6542 6.5634 6.4903 6.0820 6.2940 6.0195 5.5499 3.5931 3.6282 4.2579 4.9637 5.2469 1.0154 2.1652 4.3815 1.9166
× × × × × × × × × × × × × × × × × × × × ×
10-6 10-7 10-7 10-7 10-7 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-8 10-7 10-7 10-7 10-6
Pt error 4.0985 3.9133 3.8866 3.7391 3.3257 3.2196 2.7840 2.5691 2.1730 1.4920 1.4430 9.0957 7.0482 7.3372 1.0229 1.1961 1.8687 3.1711 3.3606 3.6753 4.5763
× × × × × × × × × × × × × × × × × × × × ×
105 105 105 105 105 105 105 105 105 105 105 104 104 104 105 105 105 105 105 105 105
8 whereas Table 4 presents the least-squares errors from experimental data for these simulations. From Figure 8 and Table 4, it can be seen that the highest deviation from experimental values is found when g ) (1 (i.e., narrow forwardly directed peak and narrow backwardly directed peak scattering). For g values in the range -0.7 < g < 0.7, the differences from MC simulations and experimental values are not very large, less than 10% in all cases. These findings suggest that for the mentioned range of asymmetry factors, a precise evaluation of the mode of reflection of the scattered photons is not very critical for a good representation of the experimental values. Results reported in this study are in agreement with those found by Pasquali et al.11 These authors studied two different distribution density functions: isotropic and diffuse phase functions. These authors concluded that both phase functions render good modeling of the radiation field in photoreactors. In the present study, the most accurate representation for the LVREA and Pt profiles were obtained with g ) 0 (isotropic scattering) for DP25 and with g ) -0.2 for Aldrich. Thus, the isotropic phase function can be used in MC simulations for both, DP25 and Aldrich catalysts. However, for Aldrich a weak backward scattering mode produced better simulation results when compared to the experimental values. These findings differ from those reported by Satuf et al.,25 who found that, for DP25, g varies from 0.6 to 0.4 and that, for Aldrich, g varies from 0.8 to 0.4 for a wavelength range 295-405 nm.
In photocatalytic systems, the slurry system contains a countless number of irregular Titania particles. However since TiO2 aggregates, this creates smoother aggregate shapes and explains the good results obtained in MC simulation with isotropic scattering averaging individual particle shapes.28 Figure 9 shows the radial profiles for the LVREA at different photocatalyst concentrations for DP25 and Aldrich for the isotropic scattering mode. It can be observed that the LVREA exhibits a quick uniform drop with the radial coordinate. One can also notice that in cases where the photocatalyst concentration is high the particles closer to the radiation source absorb most of the radiation entering the reactor. Finally, one can also conclude that if a sufficiently high photocatalyst concentration is used, a close-to-the-wall highly irradiated zone with dark areas toward the reactor center line develops. Thus, there should be an “optimum photocatalyst concentration” which maximizes conversion in the photocatalytic reactor. This optimum photocatalyst concentration also provides an optimally irradiated condition inside the reactor. Photocatalyst concentrations above this maximum show an essentially negligible effect on LVREA. According to Salaices et al.,17 and in agreement with the MC simulations, this optimum is achieved when Pt ) 0.2%Pi in a Photo-CREC Water-II unit. Furthermore, one can also notice that a great advantage of the MC simulations is that this optimum photocatalyst concentration can also be predicted by determining the LVREA at different catalyst concentrations. Additionally, one can also calculate the apparent optical thickness (τapp) which is directly related to the optimal catalyst loading.15 The apparent optical thickness for DP25 is calculated by following the procedure presented by Toepfer et al.,23 it is found that τapp ) 2.8476 for a CDP25 ) 0.14 g L-1, an annular thickness of 2.64 cm (see Figure 2), and a scattering albedo ω ) 0.89. The calculated apparent optical thickness is in agreement with the optimum value presented by Li Puma et al.15 for a plug flow reactor by using the six flux method. This comparison gives further credit to our optimal catalyst concentration. One expects that this optimum photocatalyst concentration be close to the one as that determined in photocatalytic reaction experiments. Figure 10a reports the LVREA profiles for different concentrations of DP25, the arrow shows that the LVREA for Pt ) 0.2%Pi is reached at 0.14 g L-1 of DP25. It is interesting to observe that in both cases the optimum photocatalyst concentrations are in agreement with the ones inferred from the overall reaction rate for phenol degradation, as shown in Figure 10b. It can be seen that the overall reaction rate reaches a maximum value of about 7.0 µmol-C L-1 min-1 for catalyst concentrations higher than 0.14 g L-1. As a result one can conclude that the MC simulations are not only valuable to define the optimum photocatalyst concentration leading to optimum irradiation but also an
Figure 9. Radial profiles of LVREA for different concentrations of (a) DP25 and (b) Aldrich obtained by MC: (O) 0.14, (2) 0.07, and (]) 0.04 g L-1.
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
10533
Figure 10. (a) LVREA inside the Photo-CREC Water-II as a function of DP25 loading. (b) Overall reaction rate for phenol degradation versus DP25 concentration (as presented by Salaices et al.17).
excellent tool to identify the operating conditions leading to a best possible photocatalytic rates for a given photocatalytic reactor configuration.
while developing phenol photocatalytic degradation rates experiments and with the apparent optical thickness reported in the literature.
6. Conclusions
Acknowledgment
On the basis of the results of the present study it is possible to conclude the following: (a) Monte Carlo based method can be employed to simulate the UV radiation field in an annular heterogeneous reactor using two different TiO2 photocatalysts (Aldrich and DP25). The MC method is an effective tool for solving the RTE, providing an easy to use and easy to apply alternative to circumvent the problems associated with analytical solutions. The MC simulations can be applied for virtually any reactor configuration or geometry allowing precise predictions of optimum catalyst concentrations and reactor designs. (b) MC simulations with backscattered photons reaching the BL lamp and reflected with different angular and longitudinal angles provides a more accurate result for the LVREA than when the lamp absorbs such photons. It is demonstrated that for isotropic scattering, the spectral distribution for the absorption and scattering coefficients used in the MC simulations provides a more accurate evaluation of the LVREA than when wavelength-averaged absorption and scattering coefficients are used. (c) Narrow backward and forward peaks (g ) -1 and 1, respectively) in the HG phase function are not suitable for MC simulations. On the other hand, it is demonstrated that a g value close to zero provides good representation for the experimental LVREA. It is found that for the range -0.7 < g < 0.7, differences from MC simulations and experimental values are not very large, less than 10% in all cases. This suggests that the adoption of a specific phase function is not crucial for a good representation of the radiation field, provided it is kept in the -0.7 < g < 0.7 range. It is shown comparing the light absorption rates and the transmitted radiation from both experimental observations and the MC simulations that there is satisfactory agreement. Therefore, one can use the MC simulation as an effective tool in finding the LVREA for concentric photocatalytic reactors designed on the same principles than the Photo-CREC Water-II. (d) The LVREA reaches a maximum value for DP25 concentration at the optimum photocatalyst concentration of 0.14 g L-1. It is further demonstrated that this optimum catalyst concentration for reactor irradiation is in close agreement with the optimum value found
The authors acknowledge the scholarship granted to Jesus Moreira by the Consejo Nacional de Ciencia y Tecnologı´a (CONACYT), Me´xico. We also express our appreciation to the Natural Sciences and Engineering Research Council of Canada for the financial support of this investigation. B.S. would like to especially acknowledge PIFI 2007-33-07 from Mexico and the grant CONACYT-Ciencia-Basica 2007-83144. Nomenclature Ebg ) energy band gap (eV) ex, ey, ez ) direction cosines g ) asymmetry factor G ) total incident radiation (einstein m-2 s-1) I ) radiation intensity (einstein m-2 s-1 sr-1) J ) emission term in the RTE l ) length of the photon flight (m) L ) length of the lamp (m) Pi|Cf0+ ) rate of photons transmitted when the catalyst concentration approaches to zero (einstein s-1) P(a) ) absorption criterion P(l) ) probability that the photon travels a distance l p(Ω - Ω′) ) phase function for scattering in RTE Pa ) rate of absorbed photons (einstein s-1) Pa-wall ) rate of photons absorbed by the inner Pyrex glass (einstein s-1) Pbs ) rate of backscattered photons exiting the system (einstein s-1) Pi ) rate of photons reaching the reactor inner surface (einstein s-1) Pns ) rate of transmitted nonscattered radiation (einstein s-1) Po ) rate of photons emitted by the BL lamp (einstein s-1) Pt ) rate of transmitted photons (einstein s-1) Pfs ) rate of forward-scattering radiation (einstein s-1) q(θ, z) ) net radiative flux over the lamp emission spectrum r ) radial coordinate R1-6 ) random number s ) linear coordinate along the direction Ω, (m) Wcat ) catalyst loading (g m-3) x, y, z ) Cartesian coordinates (m) Acronyms BL lamp ) black light lamp DO ) discrete ordinate
10534
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
DP25 ) TiO2 Degussa P25 FV ) finite volume LVREA ) local volumetric rate of energy absorption MC ) Monte Carlo method PTEF ) photochemical thermodynamic efficiency factor RTE ) radiative transfer equation Greek Letters β ) extinction coefficient (m-1) φ ) azimuth angle (rad), also scattering equatorial angle (rad) κ ) absorption coefficient (m-1) λ ) wavelength of radiation (m) θ ) zenith angle (rad), also scattering angular angle (rad) σ ) scattering coefficient (m-1) τapp ) apparent optical thickness ω ) scattering albedo (σ/(σ + κ)) Ω ) solid angle (steradian) Superscripts E ) denotes emission of light Subscripts λ ) indicates dependence on wavelength
Literature Cited (1) Fujishima, A.; Zhang, X.; Tryk, D. A. TiO2 photocatalysis and related surface phenomena. Surf. Sci. Rep. 2008, 63, 515. (2) Diebold, U. The surface of titanium dioxide. Surf. Sci. Rep. 2003, 48, 53. (3) de Lasa, H.; Serrano, B.; Salaices, M. Photocatalytic Reaction Engineering; Springer: New York, 2005. (4) Herrmann, J. M. Heterogeneous Photocatalysis: fundamentals and applications to the removal of various types of aqueous pollutants. Catal. Today. 1999, 53, 115. (5) Cassano, A.; Alfano, O. Reaction engineering of suspended solid heterogeneous photocatalytic reactors. Catal. Today. 2000, 58, 167. (6) Ray, A. K.; Chen, D.; Li, F. Effect of Mass Transfer and Catalyst Layer Thickness on Photocatalytic Reaction. AIChE J. 2000, 46, 1034. (7) Salaices, M.; Serrano, B.; de Lasa, H. Experimental Evaluation of photon absorption in an aqueous TiO2 Slurry reactor. Chem. Eng. J. 2002, 90, 219. (8) Sun, L.; Bolton, J. R. Determination of quantum yield for the photochemical generation of hydroxyl radicals in TiO2 suspensions. J. Phys. Chem. 1996, 100, 4127. (9) Serrano, B.; de Lasa, H. Photocatalytic Degradation of Water Organic Pollutants. Kinetic Modeling and Energy Efficiency. Ind. Eng. Chem. Res. 1997, 36, 4705. (10) Pareek, V.; Ching, S.; Tade, M.; Adesina, A. A. Light intensity distribution in heterogeneous photocatalytic reactors. Asia-Pac. J. Chem. Eng. 2008, 3, 171. (11) Pasquali, M.; Santarelli, F.; Porter, J. F.; Yue, P. L. Radiative Transfer in Photocatalytic Systems. AIChE J. 1996, 42, 532.
(12) Yokota, T.; Cesur, S.; Suzuki, H.; Baba, H.; Takahata, Y. Anisotropic Scattering Model for the Estimation of Light Absorption Rates in Photoreactor with Heterogeneous Medium. J. Chem. Eng. Jpn. 1999, 32, 314. (13) Colina-Marquez, J.; Machucha-Martinez, F.; Li Puma, G. Photocatalytic Mineralization of Commercial Herbicides in a Pilot-Scale Solar CPC Reactor: Photoreactor Modeling and Reaction Kinetics Constants Independent of Radiation Field. EnViron. Sci. Technol. 2009, 43, 8953. (14) Brucato, A.; Cassano, A. E.; Grisafi, F.; Montante, G.; Rizzuti, L.; Vella, G. Estimating Radiant Fields in Flat Heterogeneous Photoreactors by the Six-Flux Model. AIChE J. 2006, 52 (11), 3882. (15) Li Puma, G.; Brucato, A. Dimensionless analysis of slurry photocatalytic reactors using two-flux and six-flux radiation absorptionscattering models. Catal. Today. 2007, 122, 78. (16) Pareek, V. K.; Cox, S.; Adesina, A. A. Light intensity Distribution in Photocatalytic Reactors Using Finite Vol. Method. Third International Conference on CFD in the Minerals and Process Industries, CSIRO, Melbourne, Australia, 2003; p 229. (17) Salaices, M.; Serrano, B.; de Lasa, H. Photocatalytic Conversion of Organic Pollutants Extinction Coefficients and Quantum Efficiencies. Ind. Eng. Chem. Res. 2001, 40, 5455. (18) Changrani, R.; Raupp, G. B. Monte Carlo Simulation of the Radiation Field in a Reticulated Foam Photocatalytic Reactor. AIChE J. 1999, 45, 1085. (19) Brandi, R. J.; Citroni, M. A.; Alfano, O. M.; Cassano, A. E. Absolute quantum yields in photocatalytic slurry reactors. Chem. Eng. Sci. 2003, 58, 979. (20) Martin, C. A.; Baltanas, M. A.; Cassano, A. E. Photocatalytic reactors II. Quantum efficiencies allowing for scattering effects. An experimental approximation. J. Photochem. Photobiol. A 1996, 94, 173. (21) Romero, R. L.; Alfano, O. M.; Cassano, A. E. Cylindrical Photocatalytic Reactors. Radiation Absorption and Scattering Effects Produced by Suspended Fine Particles in an Annular Space. Ind. Eng. Chem. Res. 1997, 36, 3094. (22) Romero, R. L.; Alfano, O. M.; Cassano, A. E. Radiation Field in an Annular, Slurry Photocatalytic Reactor. 2. Model and Experiments. Ind. Eng. Chem. Res. 2003, 42, 2479. (23) Toepfer, B.; Gora, A.; Li Puma, G. Photocatalytic Oxidation of multicomponents solutions of herbicides: Reaction kinetics analysis with explicit photon absorption effects. Appl. Catal., B 2006, 68, 171. (24) Binzoni, T.; Leung, T. S.; Gandijbakhche, A. H.; Rufenacht, D. The Use of the Henyey-Greenstein phase function in Monte Carlo simulation in biomedical optics. Phys. Med. Biol. 2006, 51, N313. (25) Satuf, M. L.; Brandi, R. J.; Cassano, A. E.; Alfano, O. M. Experimental Method to evaluate the Optical Properties of Aqueous Titanium Dioxide Suspensions. Ind. Eng. Chem. Res. 2005, 44, 6643. (26) Marugan, J.; Van Grieken, R.; Alfano, O. M.; Cassano, A. E. Optical and Physicochemical Properties of Silica-Supported TiO2 Photocatalysts. AIChE J. 2006, 52, 832. (27) Ortiz-Gomez, A.; Serrano-Rosales, B.; Salaices, M.; de Lasa, H. Photocatalytic Oxidation of Phenol: Reaction Network, Kinetic Modeling, and Parameter Estimation. Ind. Eng. Chem. Res. 2007, 46, 7394. (28) Modest, M. F. RadiatiVe Heat Transfer, second ed.; Academic Press: California, 2003.
ReceiVed for reView February 18, 2010 ReVised manuscript receiVed April 22, 2010 Accepted April 26, 2010 IE100374F