- Email: [email protected]
A detailed kinetic model for combustion synthesis of titania from TiCl4
Combustion and Flame 156 (2009) 1764–1770
Contents lists available at ScienceDirect
Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e
A detailed kinetic model for combustion synthesis of titania from TiCl4 Richard H. West a, Raphael A. Shirley a, Markus Kraft a,*, C. Franklin Goldsmith b, William H. Green b a b
Department of Chemical Engineering and Biotechnology, University of Cambridge, New Museums Site, Cambridge, UK Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA
a r t i c l e
i n f o
Article history: Received 4 December 2008 Received in revised form 31 March 2009 Accepted 9 April 2009 Available online 6 June 2009 Keywords: Titanium dioxide Kinetic model Chemical mechanism Rapid compression machine Titanium tetrachloride
a b s t r a c t The combustion of TiCl4 to synthesize TiO2 nanoparticles is a multimillion tonne per year industrial process, the fundamental details of which are still not known. The gas-phase kinetic model presented by West et al. [R.H. West, M.S. Celnik, O.R. Inderwildi, M. Kraft, G.J.O. Beran, W.H. Green, Ind. Eng. Chem. Res. 46 (19) (2007) 6147–6156] is improved upon using density functional theory (DFT) and variational transition state theory (VTST) calculations. The pressure-dependent rate expression for the reaction TiCl3 + O2 TiO2Cl3 is found using VTST, a stable Ti2 O2 Cl6 species is located on the minimum energy pathway for TiCl3 þ TiO2 Cl3 2TiOCl3 , and a number of new elementary reactions are added. Thermochemical data are provided for Ti2 O2 Cl6 , Ti2 O2 Cl5 and TiCl2 OCl. The new kinetic model is used to simulate a rapid compression machine (RCM) and a plug flow reactor (PFR) described in the literature. Agreement with the RCM measurements is good, but simulations of the PFR are less satisfying, suggesting that surface deposition on the reactor walls may have dominated these measurements, which have been the basis of many theoretical models. Finally, the gas-phase kinetic model is coupled to a particle population balance model (PBM) incorporating inception, coagulation, growth, and sintering. Ó 2009 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
1. Introduction Titanium dioxide (TiO2) is widely used as a pigment, as a catalyst support, and as a photocatalyst. The combustion of titanium tetrachloride (TiCl4) to synthesize TiO2 nanoparticles is a multimillion tonne per year industrial process [1]. In this ‘‘chloride” process, purified titanium tetrachloride is oxidised at high temperatures (1500–2000 K) in a pure oxygen plasma or flame to produce TiO2 particles [2,3]. This paper builds on the work detailed in [4], which provides a full introduction to this field. There have been a few experimental investigations of the kinetics of the system; in the most frequently cited study, Pratsinis et al. [5] reported the overall oxidation kinetics of TiCl4 at 973–1273 K determined by measuring the concentration of unreacted TiCl4 leaving a plug flow reactor (PFR). The rate was found to be first order with respect to TiCl4 and largely independent of O2 concentration, with a rate coefficient
k ¼ 8:26 104 s1 exp
88:8kJ=mol RT
ð1Þ
where T is the temperature and R is the molar gas constant. This experimental work has been the basis for many theoretical studies of TiO2 nano-particle dynamics [6–9]. However, these * Corresponding author. E-mail address: [email protected] (M. Kraft). URL: http://como.cheng.cam.ac.uk (M. Kraft).
experiments were performed at temperatures and pressures much lower than those used in the industrial process. Extrapolation of kinetic data to this degree must be treated with caution. Attempts have also been made to measure the kinetics in a rapid compression machine (RCM), allowing higher temperatures and pressures to be reached for short times [10]. Recently a framework for a comprehensive model that attempts to describe the phenomena at different scales has been developed [11]. Ab-initio and density functional theory (DFT) investigations [4] provided the thermochemical data needed for the first thermodynamically consistent detailed kinetic model of the gas-phase combustion of TiCl4. Previous models greatly simplify the chemical processes to a single step and are unable to capture the details of temperature and concentration dependencies. Additionally the effect of H2O and other additives must eventually be included. This is only possible with a detailed kinetic model. Coupled to this detailed chemistry was a new population balance model, solved using a stochastic technique, that extended previous surface-volume models to track primary particles within each agglomerate in the population. The mechanism for surface growth of TiO2 (deposition from the gas phase) is the subject of both experimental [12] and computational [13,11] studies. Although progress is being made towards a better understanding of the surface chemistry of this system, due to a lack of detailed kinetic data, current population balance models [14,7–9,11], like this work, rely on the simple first-order rate expression measured by Ghoshtagore at 673–1020 K [15]:
0010-2180/$ - see front matter Ó 2009 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.combustflame.2009.04.011
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770
dNTiO2 dt
74:8 kJ=mol ¼ A½TiCl4 4:9 103 cm s1 exp RT
ð2Þ
where N TiO2 is the amount of TiO2 deposited, A is the surface area, and [TiCl4] is the concentration of TiCl4 in the gas phase. The purpose of this paper is (i) to improve the current detailed chemical model by identifying new species and reactions and providing thermodynamic and kinetic data; (ii) to validate the mechanism by comparing it to experiments performed in an RCM; and (iii) to simulate the PFR that led to the widely-used overall reaction rate expression, investigating the possible effect of wall deposition and the influence of temperature on particle size distributions. 2. Reaction mechanism The starting point for our reaction mechanism is the kinetic model detailed in [11]. This model includes 51 reactions, most of them with only rough estimates of their rates. The potential energy surfaces (PESs) for several reactions were investigated using the DFT package DMol3. We focus here on two reactions. First
TiCl3 þ O2 TiO2 Cl3
ðR1Þ
(R19 in [11]) because a flux analysis revealed it lies on the main reaction pathway, and it is likely to be greatly slowed by falloff. Second
TiCl3 þ TiO2 Cl3 2TiOCl3 (R21 in [11]), because a sensitivity analysis of the original kinetic model, detailed in [11], found the overall reaction rate to be most sensitive to the rate of this elementary reaction. 2.1. Reaction TiCl3+O2 TiO2Cl3 Reaction (R1) is barrierless (Fig. 1) and the rate of stabilization of TiCl3 O2 is pressure dependent. Quantum-chemical calculations were repeated using the CBS-QB3 [16] compound method in Gaussian03 [17]. The CBS-QB3 energies were calculated for the reactants and product. The B3LYP/6-311G(2d,d,p) Ti–O bond length was Re ¼1.99 Å. For the transition state minimum energy path (MEP), constrained optimizations were performed at fixed Ti–O bond lengths of R ¼2.5, 3.0, 3.5, 4.0, 4.5, and 5.0 Å., followed by CBS-QB3 calculations. The pressure-dependent reaction rate was calculated using a RRKM/Master Equation code, VariFlex [18]. The method used to calculate the variational transition state theory rate constant is similar to the method described in [19], in which a rigid-rotor harmonic-oscillator approximation is made for each step in the MEP. For both the stabilized product and the six fixed points on the MEP, the ClTi–OO torsional frequency was treated as a free rotor. Additional anharmonic effects may be present –
1765
e.g. the umbrella modes in TiCl3 O2 and the MEP with frequencies below 200 cm1 – but were not included in the calculations. A single exponential down model was used for the energy transfer due to collision, with a DEdown given by Eref ðT=298 KÞ0:8 , where Eref depends on the bath gas, here 150 cm1 for Ar and 350 cm1 for O2. No literature values could be found for the Lennard–Jones (LJ) collision parameters for the TiCl3 and TiCl3O2 complexes. Instead, the LJ collision parameters for TiCl3 were taken as the average of NH3 (since it has roughly the correct shape) and CHCl3 (since it has the correct number of Cl atoms). The LJ collision parameters for the TiCl3O2 complex were taken as the average of the bath gas and the estimated values for TiCl3. For all averaging, the LJ r was calculated from the arithmetic mean of the two values, and the LJ was calculated from the geometric mean. Although the collision parameters for TiCl3O2 are rough estimates, the final model is not expected to be very sensitive to these parameters. All LJ values for the bath gases and CHCl3 were taken from [20] and for NH3 were taken from [21]. Variflex calculations were performed at pressures between 0.001 and 30 bar, and temperatures between 400 and 2000 K. The resulting kðT; pÞ matrix was fit to a standard Lindemann [22] equation with a Troe broadening parameter [23] for use in Chemkin [24] and Cantera [25]; these parameters are given in Table 1. 2.2. Reaction TiCl3+TiO2Cl3 2TiOCl3 When investigating the PES for the exothermic reaction TiCl3 þ TiO2 Cl3 2TiOCl3 , the intermediate Ti2 O2 Cl6 was found to be stable (a minimum on the PES), with an electronic energy 225 kJ/mol below the products. In the new model this global reaction was replaced with two elementary reaction steps:
TiO2 Cl3 þ TiCl3 Ti2 O2 Cl6
ðR2Þ
2TiOCl3 Ti2 O2 Cl6
ðR3Þ
Both reactions are written here in the exothermic direction; notice that R2 and the reverse of R3 add up to the global reaction. Using a linear synchronous transit search method, and constrained geometry optimizations every 0.5 Å along the reaction coordinate, with the HCTH functional in DMol3, the minimum energy pathway (MEP) for R3 was found to be barrierless. Likewise, no transition state (saddle point on the PES) could be found for R2 at this level of theory, nor any significant barrier. As a first approximation we have assumed both reactions to be barrierless in the exothermic direction, with pre-exponential factors given by the collision limit 1013 cm3 mol1 s1. Following the addition of Ti2 O2 Cl6 to the kinetic model, a number of additional reactions were added. The abstraction of one Cl by either Cl or TiCl3 radicals leads to Ti2 O2 Cl5 , another new species; further Cl abstraction leads to Ti2 O2 Cl4 , a species already in the kinetic model. These four reactions were added to the model (R4– R7).
Energy (kJ/mol)
140
2.3. Additional reactions
120 100 80 60 40 20 0 1
2
3
4
5
6
7
Ti-O separation (Å) Fig. 1. Minimum energy pathway for reaction (R1) calculated with HCTH density functional in DMol3.
In addition to the four reactions described above, some more reactions were added to the kinetic model, including reactions of the species TiCl2 OCl which had not previously been considered. See Table 1 for all the new reactions. The rate expression for R16 was taken from [26], and that for R17 was found by fitting a modified Arrhenius expression to all the data in the NIST kinetics database [27]. The remaining rate parameters were estimated as having a pre-exponential factor, A, given by the collision limit 1 1013 cm3 mol1 s1, and the activation energy, Ea , in the exothermic direction was taken to be zero.
1766
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770
Table 1 Reaction mechanism equations. No
DH298K a
Reaction
Ab 35
277 1.92510 P 1 : TiCl3+O2 TiO2Cl3 1.0601036 c P 0 : TiCl3+O2+M TiO2Cl3+M 5 ¼ 26:93 K, T ¼ 10 K, T ¼ 5219:3 K Troe parameters: a ¼ 0:1183, T
R1
n
Ea a
6:577 6:319
41.4 0
Removed from mechanism TiCl3 þ TiO2 Cl3 2TiOCl3 Replacements R2 TiO2Cl3+TiCl3 Ti2O2Cl6 R3 2 TiOCl3 Ti2O2Cl6
7
1.001013
0
0
232 225
1.001013 1.001013
0 0
0 0
Additional reactions R4 Cl2+Ti2O2Cl5 Cl+Ti2O2Cl6 R5 Cl+Ti2O2Cl5 Cl2+Ti2O2Cl4 R6 TiCl3+Ti2O2Cl5 TiCl4+Ti2O2Cl4 R7 TiCl3+Ti2O2Cl6 TiCl4+Ti2O2Cl5 R8 TiCl2OCl TiOCl2+Cl R9 TiCl2OCl+Cl TiCl3+ClO R10 TiCl2OCl+Cl TiOCl3+Cl R11 TiCl2OCl+Cl Cl2+TiOCl2 R12 TiCl3+ClOO TiCl4+O2 R13 TiCl4+O3 TiCl3+ClO+O2 R14 TiOCl3+O3 TiO2Cl3+O2 R15 TiO2Cl2+ClOO TiO2Cl3+O2 R16 ClOO+Cl Cl2+O2 R17 O3+O 2O2
110 401 546 35 2 42 164 244 363 226 277 314 219 391
1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.391014 5.471012
0 0 0 0 0 0 0 0 0 0 0 0 0 0.00322
0 0 0 0 0 0 0 0 0 226 0 0 0 17.4
a b c
Ref.
[26] [27]
kJ mol1. cm3 mol1 s1. cm6 mol2 s1.
2.4. Thermochemistry For species added to the kinetic model in this work, thermochemical data were calculated using the same method described in [4]. Molecular geometries (Fig. 2) were optimized and analytical harmonic frequencies calculated in Gaussian 03 [17] using the B97-1 density functional [28], as recommended by Boese et al. [29]. The basis set was 6-311+G(d,p), a supplemented, triple-zeta basis set that should be large enough to ensure that the basis set truncation error is comparable with, or smaller than, the inherent errors in the DFT [29]. To establish the spin multiplicity of the ground state, the lowest-energy state of different spin multiplicity was also calculated. Ti2 O2 Cl6 is a singlet in its ground state; TiCl2 OCl and Ti2 O2 Cl5 are doublets. Isogyric reaction schemes are desired to link the calculated species to reference species with well-known enthalpies of formation. The choice is especially important, as explained in [4]. The standard enthalpy of formation of Ti2 O2 Cl6 was derived from the isogyric reaction Ti2 O2 Cl6 2TiOCl2 þ Cl2 and Ti2 O2 Cl5 from Ti2 O2 Cl5 þ TiCl4 2TiOCl2 þ Cl2 þ TiCl3 . TiCl2 OCl was derived from the anisogyric reaction TiCl2 OCl TiOCl2 þ Cl so is expected to be the least accurately determined Df H298K . Heat capacities (C p ), integrated heat capacities (HðTÞ Hð0KÞ) and entropies (S) were calculated for temperatures in the range 100–3000 K using the rigid-rotator harmonic-oscillator (RRHO)
2.206
1.777
Ti
120.48o
O
2.221
TiCl2OCl
1.665
120.30o
approximation, taking unscaled vibrational frequencies and rotational constants from the B97-1 calculations (see Table 2). Polynomials in the NASA form [30] for use in Chemkin and Cantera were fitted to C p ðTÞ=R, H and S over the temperature ranges 100 K–T c and T c –3000 K, constrained to ensure C p ðTÞ and its first two derivatives matched at the common temperature T c , which was varied to optimise the overall fit for each species. 3. Population balance modelling To simulate the comprehensive model of chemistry and particle dynamics, we use the extended surface-volume model with primary particle tracking described in [32,31,11], coupled to the gas-phase chemistry simulation using operator splitting [36]. Numerical validation of this population balance solver is provided in [36,33]. An agglomerate particle in the population balance is described by the number of TiO2 monomers in the particle, M, the surface area of the particle, A, and the number of TiO2 monomers in each of the primary particles that make up the agglomerate, (m1 ; m2 ; m3 ; . . .). All other particle properties, and hence particle process rates, are calculated from M and A. In the current study, like in [11], the collision of any two molecules containing two or more Ti atoms each is treated as a particle inception. The rate of these events is dependent on the concentration of the species involved, and is estimated by the population
118.21o 2.005 1.987
2.195 103.72o
O Ti 107.85o 43.62o O 2.172 2.209 94.61o o Ti
120.26
Ti2O2Cl5
93.76o 2.192 1.952 o 104.62 2.185 110.85o
Ti
2.175 2.047 100.36o
O O
97.05o 2.185 o 104.64 2.192
Ti
105.98o
115.80o 2.175
Ti2O2Cl6
Fig. 2. Molecular structures of TiCl2 OCl, Ti2 O2 Cl5 and Ti2 O2 Cl6 . Bond lengths are in Ångströms, and unlabelled atoms are chlorine.
1767
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770 Table 2 Calculated thermochemistry at 298 K. Species
Df H298 K (kJ/ mol)
S298 K (J/ mol K)
Rot. const. (GHz)
Vibrational frequencies (cm1)
TiCl2OCl Ti2 O2 Cl6
475 1503
396.9 562.3
1.9445, 0.9729, 0.6484 0.5965, 0.2506, 0.2497
Ti2 O2 Cl5
1272
537.5
0.7562, 0.2939, 0.2611
17.1, 44.8, 93.7, 131, 174, 341, 474, 483, 950 17.9, 18.7, 36.3, 84.8, 93.6, 101, 115, 135, 135, 137, 155, 163, 206, 299, 324, 388, 404, 455, 473, 487, 505, 509, 574, 902 14.1, 28.3, 45.8, 79.5, 82.3, 100, 118, 130, 140, 162, 194, 318, 335, 366, 401, 473, 483, 493, 501, 575, 885
balance code according to collision theory. The molecules’ collision diameters were estimated from the DFT results to be 0.65 nm (the diameter of the sphere with a volume equivalent to that occupied by the Connolly surface of Ti2O2Cl4). The rate of particle growth due to surface reactions with the gas phase was calculated using the first-order expression from [15] (Eq. (2)). Because the concentration of intermediate Ti-containing species remains low, replacing the concentration of TiCl4 in Eq. (2) with the concentration of all Ti-cointaning species would only alter the reaction rate by at most 0.5%. Sintering is modelled following the approach of [34], by assuming that the excess agglomerate surface area, over that of a spherical particle with the same mass, decays exponentially with a characteristic time taken from [35]. In order to generate shape information, to aid visualisation of the particles, the positions of the primary particles within a particle are generated in a post-processing step. 4. Simulations 4.1. Rapid compression machine The PhD thesis by Raghavan describes investigations of the kinetics of TiCl4 oxidation using a rapid compression machine (RCM) [10]. A mixture of 0.3 mol% TiCl4 and 1 mol% O2 in Argon was loaded into a cylinder at reduced pressure, then rapidly compressed causing temperature and pressure of the gas to rise quickly. The temperature was calculated from the pressure, volume, and starting temperature, according to a hot-core model in which 20% of the volume (near the walls) is compressed isothermally. The piston was designed to bounce, resulting in a short peak in temperature, found to be well approximated by a Lorentzian. The reaction time reported was the half-width at half-maximum (HWHM) of a Lorentzian fitted to the top 10% of the temperature profile (Fig. 3). The concentration of TiCl4 was measured before and after compression using FTIR spectroscopy.
This RCM was simulated here with the new gas-phase kinetic model, using Cantera and Python. In our simulation, the piston was driven sinusoidally for one stroke, with a speed and compression ratio adjusted to ensure that a Lorentzian least-squares fitted to the hottest 10% of our temperature profile had a T p and HWHM equal to Raghavan’s reported values. This was repeated for the sixteen combinations of T p and t c reported in [10]. The shape of the temperature profile from this sinusoidal compression is similar to the temperature profiles of the RCM given in [10], as shown in Fig. 3. The simulation was spatially homogeneous, so following the hot-core description in [10], the simulated conversion was scaled by 0.8 (20% of the gas remained unreacted). Because the gas is heated rapidly and for a short time the walls do not heat up, so surface reactions were neglected. Fig. 4 shows the fraction of TiCl4 consumed (1 N=N 0 ) plotted against the peak temperature reached in the reactor (T p ) from Raghavan’s thesis [10] with our simulation results superimposed. The design of the RCM meant that the peak temperature T p and the reaction time tc could not be controlled independently (a higher T p usually meant a shorter tc ). The extent of reaction obviously depends on the reaction time, so simply plotting ð1 N=N 0 Þ as a function of T p without mentioning tc is misleading. The bars on the simulation results in Fig. 4 indicate the sensitivity with respect to compression time tc : the lower limit is the result of tc ¼10 ms, the upper limit is 40 ms, and the central mark is the t c that was reported for the corresponding N=N 0 measurement in [10]. Also shown (dashed line) are the results of simulations using the onestep Arrhenius rate expression (1) from [5]. The kinetic model agrees reasonably well with the experiments without any fitting of parameters. The fraction of TiCl4 consumed during compression to a given temperature is somewhat lower in
1 Simulation 0.8
Pratsinis
1600
0.6
1−(N/N 0 )
Lorentzian
1400
Simulation
1200
Raghavan
1000
T (K)
Raghavan
0.4
0.2
800 600
0
400 200
−0.2
0
500
0
0.1
0.2
0.3
time (ms) Fig. 3. Comparison of temperature profiles from the simulation and from Raghavan [10]. Included is Raghavan’s fitted lorentzian from which peak temperatures and residence times are calculated.
1000
1500
2000
T p (K) Fig. 4. Fraction of TiCl4 consumed after rapid compression in the RCM, plotted against the peak temperature reached. The lower limit shows a simulation with tc ¼10 ms, the upper limit 40 ms, and the filled diamond the tc that was reported for the corresponding measurement (white circle) from [10].
1768
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770
our simulations than the reported measurements, but, although no error bars were reported for the measurements, we expect the uncertainties in the temperature estimations to be significant (see Section 5).
higher than that observed in [5]. This may be partially because mass transport limitations were not modelled. Also, Ghoshtaghore’s rate expression is for much higher O2 partial pressures, so would be expected to over-predict wall reactions in the conditions modelled here.
4.2. Plug Flow Reactor (PFR) 4.3. Particle properties The combined simulation coupling the gas-phase kinetic model to the particle dynamics model was used to simulate the tubular flow reactor used by Pratsinis et al. [5] to measure the overall reaction kinetics. In [5], a flow of 0.2 mol% TiCl4 and 1 mol% O2 in Argon was fed to a 3.18 mm-I.D. tube, a 30 cm section of which was heated to 973–1273 K in a furnace. The concentration of TiCl4 leaving the flow tube (C o ) was measured using FTIR spectroscopy. Assuming that the reaction rate is first order in [TiCl4] with Arrhenius kinetics (keff ¼ A expðEa =RTÞ) and that all reaction occurs during the time t that the gas spends in an isothermal zone at temperature T then
ðln ðC o =C i ÞÞ=t ¼ keff ¼ A expðEa =RTÞ:
ð3Þ
A plot of ln ððlnðC o =C i ÞÞ=t Þ versus 1=T (Fig. 5) would thus give a straight line from which the apparent activation energy Ea and pre-exponential constant A can be determined. This is the analysis performed in [5], and our simulation results are processed the same way for easy comparison. The assumption of constant T for a time t was only made for the analysis: the simulations covered the entire temperature profile (given in [5]) and included thermal expansion of the gas as the temperature varied. The simulations using our new gas-phase kinetic model (Fig. 5, solid lines) consume much less TiCl4 than the observations reported in [5]. At 1133 K the keff fitted to the simulation is between 14 and 45 times lower (slower rate) than that fitted to the experiment, depending on which residence time is simulated. The apparent activation energy is higher in the simulation, so at 993 K the expt sim discrepancy increases to keff =keff ¼ 5:2 103 in the worst case. Possible reasons for this discrepancy will be discussed in a later section. The simulations were repeated allowing for deposition and reaction of TiCl4 on the inside surface of the tube, using the expression for surface growth of a TiO2 film given by Ghoshtagore [15] (Eq. (2)). In these simulations (Fig. 5, dotted lines), 500 times as much TiCl4 was consumed through reaction on the walls as through reactions in the gas phase, and the overall consumption of TiCl4 was
The flow tube reactor was simulated using the coupled population balance model. This allows access to a wealth of information far beyond the rate of consumption of TiCl4. Fig. 6 shows the size distribution of titania particles at the end of the simulation (40 cm downstream from the midpoint of the tube, when the temperature has returned to 298 K) and Fig. 7 shows TEM-style images of some of the particles. These simulated TEMs are created from the list of primary particles that is produced by the population balance code as per the method described in [36]. As the temperature of the PFR is increased, so does the average size of the particles produced. This is the same as the trend observed by [37] in a similar PFR. 4.4. Comparison of rates Although the complex pressure-dependent chemistry of this system cannot be fully described by a first-order rate expression it can still be interesting to compare the different measurements and predictions in terms of an apparent first-order rate coefficient. The determination of an effective rate coefficient from the plug flow reactor experiments was described in Section 4.2. The best fit Arrhenius expressions reported in [5] for four different O2/TiCl4 ratios are plotted on an Arrhenius plot in Fig. 8. Comparable rate coefficients for the rapid compression machine and the detailed kinetic model were found as follows. 4.4.1. From the detailed kinetic model A batch reactor at constant temperature and pressure and initial composition of 1 mol% TiCl4 and 5 mol% O2 in Argon was simulated for a time t ¼ minðthalf ; 1sÞ, which is the lesser of the time at which C o =C i ¼ 0:5 or 1 s. The choice of maximum reaction time (1 s) has only a small effect on the results. The apparent first-order rate coefficient keff was determined according to
keff ¼ ðln ðC o =C i ÞÞ=t:
ð4Þ
The simulations were repeated at pressures from 0.1 to 100 bar to investigate the pressure-dependence of the model. As shown in
4 2
100 T=993 K T=1073 K 10
−2 d P/ d log Dp
ln(k eff )
0
−4 −6
0.49 s
−8
0.75 s
T=1133 K T=1293 K
1
0.1
1.1 s −10 0.85
0.9
0.95
1.0
103 K/T Fig. 5. Arrhenius plot of plug flow reactor at three different residence times. Unconnected points are measurements reported in [5], lines are our simulations: dotted with deposition on the reactor walls and solid without.
0.01 0.0001
0.001
0.01 Dp (µm)
0.1
1
Fig. 6. Particle size distributions from the PFR simulation with residence time in the hot-zone of 0.49 s.
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770
1769
Fig. 7. Simulated TEM Images from the exit of the PFR simulation with residence time in the hot-zone of 0.49 s and furnace temperatures 993, 1017, 1073 and 1293 K (left to right).
Fig. 8, the apparent rate constant is much higher at higher pressures, increasing almost linearly with pressure. 4.4.2. From Raghavan’s RCM data To determine an effective first-order rate coefficient from the RCM measurements requires some assumptions. First, that the temperature follows a Lorentzian with HMHM of t c and peak temperature of T p . Second, that the temperature-dependence of the rate coefficient follows expðEa =RTÞ. Ea was taken to be 108 kJ/ mol as given by Raghavan, but the results are not very sensitive to this choice. With these assumptions the effective rate coefficient was determined from a single compression experiment according to
keff ¼
! ln ððN=N0 f Þ=ð1 f ÞÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : tc pRT p =Ea
ð5Þ
Where f is the fraction of the gas that is compressed isothermally according to the hot-core model. There are some difficulties with determining keff from the observed ðN=N 0 Þ values: first, if N=N0 P 1 then the apparent rate is negative; second, if N=N 0 6 f then the apparent rate is greater than infinity. The uncertainties of keff values from experiments with high and low conversions are therefore very large, making lnðkeff Þ positive and negative infinity, respectively. To illustrate this the bars in Fig. 8 shows the error in lnðkeff Þ that would result from an error of 0.1 in the measurement of ðN=N0 Þ. 5. Discussion A rapid compression machine seems to be a useful tool for measuring gas-phase kinetics at high temperatures and pressures, but at less expense than a shock-tube. However, uncertainties in the temperature measurements are likely to be significant. The peak temperature T p was deduced from measurements of pressure, volume, and the starting temperature T 0 , giving three causes for con-
Our Simulation 10
Raghavan RCM Pratsinis PFR
ln(k)
5 0 -5
100 bar 10 bar 1 bar 0.1 bar
-10 0.5
1.0
1.5
cern. First, an error in measuring T0 is amplified in T p . Second, the deduced T p is very sensitive to errors in measuring the maximum piston displacement. Third, the assumptions necessary for this deduction are hard to verify. Unfortunately, Raghavan does not give estimates of the errors in N=N 0 or T p in Fig. 4, but we expect they are not small enough to justify adjusting parameters in our kinetic model to fit these experiments. The 20% dead-zone in the model ensures that no simulation will ever match the observations in which ð1 N=N 0 Þ > 0:8, whichever kinetics are used. Plug flow reactors are often used to measure rates, and the careful measurements in [5] seem like a useful check for any kinetic model of TiCl4 oxidation. At first glance the disagreement of four orders of magnitude between the first-order rate fitted to the experiments and that fitted to the simulations suggests a problem with the new gas-phase kinetic model used in the simulations. Indeed, this may be the case, and further work on the model is recommended, especially the initiating reactions. However, there is an alternative possibility. Due to the small internal diameter of the tube, the surface-area to volume ratio is high, and including surface reaction on the tube wall has a large effect on the simulations. It is possible that some of the observed consumption of TiCl4 occurred on the walls of the tube, rather than in the gas phase. The rate expression derived from these measurements, however, has been widely used to simulate gas-phase reactions, and has often been used as an inception rate for TiO2 particle population balance models. Yoon et al. [38] observed that, in their tubular reactor of diameter 2.7 cm, up to 92% of TiO2 production was on the walls; although a different experiment, this illustrates that in at least some conditions, wall reactions can dominate. Further investigation of both gas and surface kinetics may be required to resolve this issue. Experimental investigations into the influence of tube diameter on the overall rate of decomposition could be particularly instructive. The uncertainty in the kinetic model is difficult to assess fully and a proper treatment is beyond the scope of this paper. Spectral expansion techniques such as those in Sheen et al. [39], Russi et al. [40] and Xiu et al. [41] can be used to rigorously quantify uncertainties in the kinetic rate parameters. Further work using techniques like this will enable the determination of the overall errors associated with the model. As shown in Fig. 8 the detailed kinetic model predicts that the overall reaction rate is pressure-dependent, and the temperaturedependence is too complicated to describe with a single Arrhenius expression. This decreases the accuracy of rate expressions determined experimentally over a narrow range of temperature and pressure, especially when extrapolated beyond that range, for example to model industrial reactor conditions. The current detailed model is much more versatile in terms of operating conditions. 6. Conclusions
(1000 K) / T Fig. 8. Arrhenius plots of our model at different pressures with Pratsinis [5] and Raghavan [10] for comparison.
The first detailed gas-phase kinetic model for TiCl4 oxidation [11] has been improved using DFT and VTST calculations, which
1770
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770
have provided new thermochemical and kinetic data. This improved kinetic model is in agreement with the experimental data from a rapid compression machine [10], with the results of the model falling within the estimated uncertainties of the experimental measurements. Simulations of a plug flow reactor, however, suggest that surface reactions on the reactor walls may have some impact on the experimental measurements [5]. Many models simulating TiO2 formation from the gas phase have relied on rate expressions derived from these experimental measurements. The improved kinetic model thus raises an important question in the field of modelling TiO2 particle formation. Finally, the improved gas-phase kinetic model was coupled to a particle population balance model (PBM) to investigate the influence of PFR temperature on particle size distributions. The traditional surface-volume PBM was extended to track the growth and sintering of primary particles within each particle. This allowed creation of TEM-style images of the agglomerate particles from the simulation, which are presented herein. Acknowledgments The authors would like to thank Prof. Stephen Klippenstein for help with the Variflex calculations, Tioxide Europe Limited (TEL) for the financial support of R.A.S., the EPSRC for Grant No. EP/ E01724X/1, and both the EPSRC and TEL for the financial support of R.H.W. through an Industrial CASE studentship. Appendix A. Supplementary data
[13] [14] [15] [16] [17]
[18]
[19] [20] [21] [22] [23] [24]
[25]
[26]
Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.combustflame.2009.04.011. [27]
References [1] J. Emsley, Molecules at an Exhibition, Oxford University Press, 1999. [2] R.A. Gonzalez, C.D. Musick, J.N. Tilton, Process for controlling agglomeration in the manufacture of TiO2, US Patent 5,508,015, 1996. [3] J.C. Deberry, M. Robinson, M. Pomponi, A.J. Beach, Y. Xiong, K. Akhtar, Controlled vapor phase oxidation of titanium tetrachloride to manufacture titanium dioxide, US Patent 6,387,347, 2002. [4] R.H. West, G.J.O. Beran, W.H. Green, M. Kraft, J. Phys. Chem. A 111 (18) (2007) 3560–3565. [5] S.E. Pratsinis, H. Bai, P. Biswas, M. Frenklach, S.V.R. Mastrangelo, J. Am. Ceram. Soc. 73 (1990) 2158–2162. [6] P. Koukkari, J. Niemela, Comput. Chem. Eng. 21 (1997) 245–253. [7] P.T. Spicer, O. Chaoul, S. Tsantilis, S.E. Pratsinis, J. Aerosol Sci. 33 (2002) 17–34. [8] S. Tsantilis, S.E. Pratsinis, J. Aerosol. Sci. 35 (2004) 405–420. [9] N.M. Morgan, C.G. Wells, M.J. Goodson, M. Kraft, W. Wagner, J. Comput. Phys. 211 (2) (2006) 638–658. [10] R. Raghavan, Measurement of the High-Temperature Kinetics of Titanium Tetrachloride (TiCl4) Reactions in a Rapid Compression Machine, Ph.D. thesis, Case Western Reserve University, 2001. Available from:. [11] R.H. West, M.S. Celnik, O.R. Inderwildi, M. Kraft, G.J.O. Beran, W.H. Green, Ind. Eng. Chem. Res. 46 (19) (2007) 6147–6156. [12] R.D. Smith, R.A. Bennett, M. Bowker, Phys. Rev. B: Condens. Matter Mater. Phys. 66 (3) (2002) 035409.
[28] [29] [30]
[31] [32] [33] [34] [35] [36]
[37] [38] [39] [40] [41]
O. Inderwildi, M. Kraft, ChemPhysChem 8 (2007) 444–451. S.E. Pratsinis, P.T. Spicer, Chem. Eng. Sci. 53 (1998) 1861–1868. R.N. Ghoshtagore, J. Electrochem. Soc. 117 (1970) 529–534. J.A. Montgomery Jr., M.J. Frisch, J.W. Ochterski, G.A. Petersson, J. Chem. Phys. 112 (15) (2000) 6532–6542. M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, J.A. Montgomery, Jr., T. Vreven, K.N. Kudin, J.C. Burant, J.M. Millam, S.S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G.A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J.E. Knox, H.P. Hratchian, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, P.Y. Ayala, K. Morokuma, G.A. Voth, P. Salvador, J.J. Dannenberg, V.G. Zakrzewski, S. Dapprich, A.D. Daniels, M.C. Strain, O. Farkas, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J.V. Ortiz, Q. Cui, A.G. Baboul, S. Clifford, J. Cioslowski, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, C. Gonzalez, J.A. Pople, Gaussian 03, Revision C.02, Gaussian, Inc., Wallingford, CT, 2004 (2003). S.J. Klippenstein, A.F. Wagner, R.C. Dunbar, D.M. Wardlaw, S.H. Robertson, J.A. Miller, VariFlex, software package, 2002. Available from:. W. DeMore, S. Sander, D. Golden, R. Hampson, M. Kurylo, C. Howard, A. Ravishankara, C. Kolb, M. Molina, Chemical kinetics and photochemical data for use in stratospheric modeling, Evaluation Number 15, Tech. Rep. JPL 06-2, NASA Jet Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, California 91109-8099 USA, 1997. NIST, Chemical Kinetics Database on the Web. Standard Reference Database 17, Version 7.0 (Web Version), Release 1.4, 2007. Available from:
Contents lists available at ScienceDirect
Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e
A detailed kinetic model for combustion synthesis of titania from TiCl4 Richard H. West a, Raphael A. Shirley a, Markus Kraft a,*, C. Franklin Goldsmith b, William H. Green b a b
Department of Chemical Engineering and Biotechnology, University of Cambridge, New Museums Site, Cambridge, UK Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA
a r t i c l e
i n f o
Article history: Received 4 December 2008 Received in revised form 31 March 2009 Accepted 9 April 2009 Available online 6 June 2009 Keywords: Titanium dioxide Kinetic model Chemical mechanism Rapid compression machine Titanium tetrachloride
a b s t r a c t The combustion of TiCl4 to synthesize TiO2 nanoparticles is a multimillion tonne per year industrial process, the fundamental details of which are still not known. The gas-phase kinetic model presented by West et al. [R.H. West, M.S. Celnik, O.R. Inderwildi, M. Kraft, G.J.O. Beran, W.H. Green, Ind. Eng. Chem. Res. 46 (19) (2007) 6147–6156] is improved upon using density functional theory (DFT) and variational transition state theory (VTST) calculations. The pressure-dependent rate expression for the reaction TiCl3 + O2 TiO2Cl3 is found using VTST, a stable Ti2 O2 Cl6 species is located on the minimum energy pathway for TiCl3 þ TiO2 Cl3 2TiOCl3 , and a number of new elementary reactions are added. Thermochemical data are provided for Ti2 O2 Cl6 , Ti2 O2 Cl5 and TiCl2 OCl. The new kinetic model is used to simulate a rapid compression machine (RCM) and a plug flow reactor (PFR) described in the literature. Agreement with the RCM measurements is good, but simulations of the PFR are less satisfying, suggesting that surface deposition on the reactor walls may have dominated these measurements, which have been the basis of many theoretical models. Finally, the gas-phase kinetic model is coupled to a particle population balance model (PBM) incorporating inception, coagulation, growth, and sintering. Ó 2009 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
1. Introduction Titanium dioxide (TiO2) is widely used as a pigment, as a catalyst support, and as a photocatalyst. The combustion of titanium tetrachloride (TiCl4) to synthesize TiO2 nanoparticles is a multimillion tonne per year industrial process [1]. In this ‘‘chloride” process, purified titanium tetrachloride is oxidised at high temperatures (1500–2000 K) in a pure oxygen plasma or flame to produce TiO2 particles [2,3]. This paper builds on the work detailed in [4], which provides a full introduction to this field. There have been a few experimental investigations of the kinetics of the system; in the most frequently cited study, Pratsinis et al. [5] reported the overall oxidation kinetics of TiCl4 at 973–1273 K determined by measuring the concentration of unreacted TiCl4 leaving a plug flow reactor (PFR). The rate was found to be first order with respect to TiCl4 and largely independent of O2 concentration, with a rate coefficient
k ¼ 8:26 104 s1 exp
88:8kJ=mol RT
ð1Þ
where T is the temperature and R is the molar gas constant. This experimental work has been the basis for many theoretical studies of TiO2 nano-particle dynamics [6–9]. However, these * Corresponding author. E-mail address: [email protected] (M. Kraft). URL: http://como.cheng.cam.ac.uk (M. Kraft).
experiments were performed at temperatures and pressures much lower than those used in the industrial process. Extrapolation of kinetic data to this degree must be treated with caution. Attempts have also been made to measure the kinetics in a rapid compression machine (RCM), allowing higher temperatures and pressures to be reached for short times [10]. Recently a framework for a comprehensive model that attempts to describe the phenomena at different scales has been developed [11]. Ab-initio and density functional theory (DFT) investigations [4] provided the thermochemical data needed for the first thermodynamically consistent detailed kinetic model of the gas-phase combustion of TiCl4. Previous models greatly simplify the chemical processes to a single step and are unable to capture the details of temperature and concentration dependencies. Additionally the effect of H2O and other additives must eventually be included. This is only possible with a detailed kinetic model. Coupled to this detailed chemistry was a new population balance model, solved using a stochastic technique, that extended previous surface-volume models to track primary particles within each agglomerate in the population. The mechanism for surface growth of TiO2 (deposition from the gas phase) is the subject of both experimental [12] and computational [13,11] studies. Although progress is being made towards a better understanding of the surface chemistry of this system, due to a lack of detailed kinetic data, current population balance models [14,7–9,11], like this work, rely on the simple first-order rate expression measured by Ghoshtagore at 673–1020 K [15]:
0010-2180/$ - see front matter Ó 2009 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.combustflame.2009.04.011
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770
dNTiO2 dt
74:8 kJ=mol ¼ A½TiCl4 4:9 103 cm s1 exp RT
ð2Þ
where N TiO2 is the amount of TiO2 deposited, A is the surface area, and [TiCl4] is the concentration of TiCl4 in the gas phase. The purpose of this paper is (i) to improve the current detailed chemical model by identifying new species and reactions and providing thermodynamic and kinetic data; (ii) to validate the mechanism by comparing it to experiments performed in an RCM; and (iii) to simulate the PFR that led to the widely-used overall reaction rate expression, investigating the possible effect of wall deposition and the influence of temperature on particle size distributions. 2. Reaction mechanism The starting point for our reaction mechanism is the kinetic model detailed in [11]. This model includes 51 reactions, most of them with only rough estimates of their rates. The potential energy surfaces (PESs) for several reactions were investigated using the DFT package DMol3. We focus here on two reactions. First
TiCl3 þ O2 TiO2 Cl3
ðR1Þ
(R19 in [11]) because a flux analysis revealed it lies on the main reaction pathway, and it is likely to be greatly slowed by falloff. Second
TiCl3 þ TiO2 Cl3 2TiOCl3 (R21 in [11]), because a sensitivity analysis of the original kinetic model, detailed in [11], found the overall reaction rate to be most sensitive to the rate of this elementary reaction. 2.1. Reaction TiCl3+O2 TiO2Cl3 Reaction (R1) is barrierless (Fig. 1) and the rate of stabilization of TiCl3 O2 is pressure dependent. Quantum-chemical calculations were repeated using the CBS-QB3 [16] compound method in Gaussian03 [17]. The CBS-QB3 energies were calculated for the reactants and product. The B3LYP/6-311G(2d,d,p) Ti–O bond length was Re ¼1.99 Å. For the transition state minimum energy path (MEP), constrained optimizations were performed at fixed Ti–O bond lengths of R ¼2.5, 3.0, 3.5, 4.0, 4.5, and 5.0 Å., followed by CBS-QB3 calculations. The pressure-dependent reaction rate was calculated using a RRKM/Master Equation code, VariFlex [18]. The method used to calculate the variational transition state theory rate constant is similar to the method described in [19], in which a rigid-rotor harmonic-oscillator approximation is made for each step in the MEP. For both the stabilized product and the six fixed points on the MEP, the ClTi–OO torsional frequency was treated as a free rotor. Additional anharmonic effects may be present –
1765
e.g. the umbrella modes in TiCl3 O2 and the MEP with frequencies below 200 cm1 – but were not included in the calculations. A single exponential down model was used for the energy transfer due to collision, with a DEdown given by Eref ðT=298 KÞ0:8 , where Eref depends on the bath gas, here 150 cm1 for Ar and 350 cm1 for O2. No literature values could be found for the Lennard–Jones (LJ) collision parameters for the TiCl3 and TiCl3O2 complexes. Instead, the LJ collision parameters for TiCl3 were taken as the average of NH3 (since it has roughly the correct shape) and CHCl3 (since it has the correct number of Cl atoms). The LJ collision parameters for the TiCl3O2 complex were taken as the average of the bath gas and the estimated values for TiCl3. For all averaging, the LJ r was calculated from the arithmetic mean of the two values, and the LJ was calculated from the geometric mean. Although the collision parameters for TiCl3O2 are rough estimates, the final model is not expected to be very sensitive to these parameters. All LJ values for the bath gases and CHCl3 were taken from [20] and for NH3 were taken from [21]. Variflex calculations were performed at pressures between 0.001 and 30 bar, and temperatures between 400 and 2000 K. The resulting kðT; pÞ matrix was fit to a standard Lindemann [22] equation with a Troe broadening parameter [23] for use in Chemkin [24] and Cantera [25]; these parameters are given in Table 1. 2.2. Reaction TiCl3+TiO2Cl3 2TiOCl3 When investigating the PES for the exothermic reaction TiCl3 þ TiO2 Cl3 2TiOCl3 , the intermediate Ti2 O2 Cl6 was found to be stable (a minimum on the PES), with an electronic energy 225 kJ/mol below the products. In the new model this global reaction was replaced with two elementary reaction steps:
TiO2 Cl3 þ TiCl3 Ti2 O2 Cl6
ðR2Þ
2TiOCl3 Ti2 O2 Cl6
ðR3Þ
Both reactions are written here in the exothermic direction; notice that R2 and the reverse of R3 add up to the global reaction. Using a linear synchronous transit search method, and constrained geometry optimizations every 0.5 Å along the reaction coordinate, with the HCTH functional in DMol3, the minimum energy pathway (MEP) for R3 was found to be barrierless. Likewise, no transition state (saddle point on the PES) could be found for R2 at this level of theory, nor any significant barrier. As a first approximation we have assumed both reactions to be barrierless in the exothermic direction, with pre-exponential factors given by the collision limit 1013 cm3 mol1 s1. Following the addition of Ti2 O2 Cl6 to the kinetic model, a number of additional reactions were added. The abstraction of one Cl by either Cl or TiCl3 radicals leads to Ti2 O2 Cl5 , another new species; further Cl abstraction leads to Ti2 O2 Cl4 , a species already in the kinetic model. These four reactions were added to the model (R4– R7).
Energy (kJ/mol)
140
2.3. Additional reactions
120 100 80 60 40 20 0 1
2
3
4
5
6
7
Ti-O separation (Å) Fig. 1. Minimum energy pathway for reaction (R1) calculated with HCTH density functional in DMol3.
In addition to the four reactions described above, some more reactions were added to the kinetic model, including reactions of the species TiCl2 OCl which had not previously been considered. See Table 1 for all the new reactions. The rate expression for R16 was taken from [26], and that for R17 was found by fitting a modified Arrhenius expression to all the data in the NIST kinetics database [27]. The remaining rate parameters were estimated as having a pre-exponential factor, A, given by the collision limit 1 1013 cm3 mol1 s1, and the activation energy, Ea , in the exothermic direction was taken to be zero.
1766
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770
Table 1 Reaction mechanism equations. No
DH298K a
Reaction
Ab 35
277 1.92510 P 1 : TiCl3+O2 TiO2Cl3 1.0601036 c P 0 : TiCl3+O2+M TiO2Cl3+M 5 ¼ 26:93 K, T ¼ 10 K, T ¼ 5219:3 K Troe parameters: a ¼ 0:1183, T
R1
n
Ea a
6:577 6:319
41.4 0
Removed from mechanism TiCl3 þ TiO2 Cl3 2TiOCl3 Replacements R2 TiO2Cl3+TiCl3 Ti2O2Cl6 R3 2 TiOCl3 Ti2O2Cl6
7
1.001013
0
0
232 225
1.001013 1.001013
0 0
0 0
Additional reactions R4 Cl2+Ti2O2Cl5 Cl+Ti2O2Cl6 R5 Cl+Ti2O2Cl5 Cl2+Ti2O2Cl4 R6 TiCl3+Ti2O2Cl5 TiCl4+Ti2O2Cl4 R7 TiCl3+Ti2O2Cl6 TiCl4+Ti2O2Cl5 R8 TiCl2OCl TiOCl2+Cl R9 TiCl2OCl+Cl TiCl3+ClO R10 TiCl2OCl+Cl TiOCl3+Cl R11 TiCl2OCl+Cl Cl2+TiOCl2 R12 TiCl3+ClOO TiCl4+O2 R13 TiCl4+O3 TiCl3+ClO+O2 R14 TiOCl3+O3 TiO2Cl3+O2 R15 TiO2Cl2+ClOO TiO2Cl3+O2 R16 ClOO+Cl Cl2+O2 R17 O3+O 2O2
110 401 546 35 2 42 164 244 363 226 277 314 219 391
1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.001013 1.391014 5.471012
0 0 0 0 0 0 0 0 0 0 0 0 0 0.00322
0 0 0 0 0 0 0 0 0 226 0 0 0 17.4
a b c
Ref.
[26] [27]
kJ mol1. cm3 mol1 s1. cm6 mol2 s1.
2.4. Thermochemistry For species added to the kinetic model in this work, thermochemical data were calculated using the same method described in [4]. Molecular geometries (Fig. 2) were optimized and analytical harmonic frequencies calculated in Gaussian 03 [17] using the B97-1 density functional [28], as recommended by Boese et al. [29]. The basis set was 6-311+G(d,p), a supplemented, triple-zeta basis set that should be large enough to ensure that the basis set truncation error is comparable with, or smaller than, the inherent errors in the DFT [29]. To establish the spin multiplicity of the ground state, the lowest-energy state of different spin multiplicity was also calculated. Ti2 O2 Cl6 is a singlet in its ground state; TiCl2 OCl and Ti2 O2 Cl5 are doublets. Isogyric reaction schemes are desired to link the calculated species to reference species with well-known enthalpies of formation. The choice is especially important, as explained in [4]. The standard enthalpy of formation of Ti2 O2 Cl6 was derived from the isogyric reaction Ti2 O2 Cl6 2TiOCl2 þ Cl2 and Ti2 O2 Cl5 from Ti2 O2 Cl5 þ TiCl4 2TiOCl2 þ Cl2 þ TiCl3 . TiCl2 OCl was derived from the anisogyric reaction TiCl2 OCl TiOCl2 þ Cl so is expected to be the least accurately determined Df H298K . Heat capacities (C p ), integrated heat capacities (HðTÞ Hð0KÞ) and entropies (S) were calculated for temperatures in the range 100–3000 K using the rigid-rotator harmonic-oscillator (RRHO)
2.206
1.777
Ti
120.48o
O
2.221
TiCl2OCl
1.665
120.30o
approximation, taking unscaled vibrational frequencies and rotational constants from the B97-1 calculations (see Table 2). Polynomials in the NASA form [30] for use in Chemkin and Cantera were fitted to C p ðTÞ=R, H and S over the temperature ranges 100 K–T c and T c –3000 K, constrained to ensure C p ðTÞ and its first two derivatives matched at the common temperature T c , which was varied to optimise the overall fit for each species. 3. Population balance modelling To simulate the comprehensive model of chemistry and particle dynamics, we use the extended surface-volume model with primary particle tracking described in [32,31,11], coupled to the gas-phase chemistry simulation using operator splitting [36]. Numerical validation of this population balance solver is provided in [36,33]. An agglomerate particle in the population balance is described by the number of TiO2 monomers in the particle, M, the surface area of the particle, A, and the number of TiO2 monomers in each of the primary particles that make up the agglomerate, (m1 ; m2 ; m3 ; . . .). All other particle properties, and hence particle process rates, are calculated from M and A. In the current study, like in [11], the collision of any two molecules containing two or more Ti atoms each is treated as a particle inception. The rate of these events is dependent on the concentration of the species involved, and is estimated by the population
118.21o 2.005 1.987
2.195 103.72o
O Ti 107.85o 43.62o O 2.172 2.209 94.61o o Ti
120.26
Ti2O2Cl5
93.76o 2.192 1.952 o 104.62 2.185 110.85o
Ti
2.175 2.047 100.36o
O O
97.05o 2.185 o 104.64 2.192
Ti
105.98o
115.80o 2.175
Ti2O2Cl6
Fig. 2. Molecular structures of TiCl2 OCl, Ti2 O2 Cl5 and Ti2 O2 Cl6 . Bond lengths are in Ångströms, and unlabelled atoms are chlorine.
1767
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770 Table 2 Calculated thermochemistry at 298 K. Species
Df H298 K (kJ/ mol)
S298 K (J/ mol K)
Rot. const. (GHz)
Vibrational frequencies (cm1)
TiCl2OCl Ti2 O2 Cl6
475 1503
396.9 562.3
1.9445, 0.9729, 0.6484 0.5965, 0.2506, 0.2497
Ti2 O2 Cl5
1272
537.5
0.7562, 0.2939, 0.2611
17.1, 44.8, 93.7, 131, 174, 341, 474, 483, 950 17.9, 18.7, 36.3, 84.8, 93.6, 101, 115, 135, 135, 137, 155, 163, 206, 299, 324, 388, 404, 455, 473, 487, 505, 509, 574, 902 14.1, 28.3, 45.8, 79.5, 82.3, 100, 118, 130, 140, 162, 194, 318, 335, 366, 401, 473, 483, 493, 501, 575, 885
balance code according to collision theory. The molecules’ collision diameters were estimated from the DFT results to be 0.65 nm (the diameter of the sphere with a volume equivalent to that occupied by the Connolly surface of Ti2O2Cl4). The rate of particle growth due to surface reactions with the gas phase was calculated using the first-order expression from [15] (Eq. (2)). Because the concentration of intermediate Ti-containing species remains low, replacing the concentration of TiCl4 in Eq. (2) with the concentration of all Ti-cointaning species would only alter the reaction rate by at most 0.5%. Sintering is modelled following the approach of [34], by assuming that the excess agglomerate surface area, over that of a spherical particle with the same mass, decays exponentially with a characteristic time taken from [35]. In order to generate shape information, to aid visualisation of the particles, the positions of the primary particles within a particle are generated in a post-processing step. 4. Simulations 4.1. Rapid compression machine The PhD thesis by Raghavan describes investigations of the kinetics of TiCl4 oxidation using a rapid compression machine (RCM) [10]. A mixture of 0.3 mol% TiCl4 and 1 mol% O2 in Argon was loaded into a cylinder at reduced pressure, then rapidly compressed causing temperature and pressure of the gas to rise quickly. The temperature was calculated from the pressure, volume, and starting temperature, according to a hot-core model in which 20% of the volume (near the walls) is compressed isothermally. The piston was designed to bounce, resulting in a short peak in temperature, found to be well approximated by a Lorentzian. The reaction time reported was the half-width at half-maximum (HWHM) of a Lorentzian fitted to the top 10% of the temperature profile (Fig. 3). The concentration of TiCl4 was measured before and after compression using FTIR spectroscopy.
This RCM was simulated here with the new gas-phase kinetic model, using Cantera and Python. In our simulation, the piston was driven sinusoidally for one stroke, with a speed and compression ratio adjusted to ensure that a Lorentzian least-squares fitted to the hottest 10% of our temperature profile had a T p and HWHM equal to Raghavan’s reported values. This was repeated for the sixteen combinations of T p and t c reported in [10]. The shape of the temperature profile from this sinusoidal compression is similar to the temperature profiles of the RCM given in [10], as shown in Fig. 3. The simulation was spatially homogeneous, so following the hot-core description in [10], the simulated conversion was scaled by 0.8 (20% of the gas remained unreacted). Because the gas is heated rapidly and for a short time the walls do not heat up, so surface reactions were neglected. Fig. 4 shows the fraction of TiCl4 consumed (1 N=N 0 ) plotted against the peak temperature reached in the reactor (T p ) from Raghavan’s thesis [10] with our simulation results superimposed. The design of the RCM meant that the peak temperature T p and the reaction time tc could not be controlled independently (a higher T p usually meant a shorter tc ). The extent of reaction obviously depends on the reaction time, so simply plotting ð1 N=N 0 Þ as a function of T p without mentioning tc is misleading. The bars on the simulation results in Fig. 4 indicate the sensitivity with respect to compression time tc : the lower limit is the result of tc ¼10 ms, the upper limit is 40 ms, and the central mark is the t c that was reported for the corresponding N=N 0 measurement in [10]. Also shown (dashed line) are the results of simulations using the onestep Arrhenius rate expression (1) from [5]. The kinetic model agrees reasonably well with the experiments without any fitting of parameters. The fraction of TiCl4 consumed during compression to a given temperature is somewhat lower in
1 Simulation 0.8
Pratsinis
1600
0.6
1−(N/N 0 )
Lorentzian
1400
Simulation
1200
Raghavan
1000
T (K)
Raghavan
0.4
0.2
800 600
0
400 200
−0.2
0
500
0
0.1
0.2
0.3
time (ms) Fig. 3. Comparison of temperature profiles from the simulation and from Raghavan [10]. Included is Raghavan’s fitted lorentzian from which peak temperatures and residence times are calculated.
1000
1500
2000
T p (K) Fig. 4. Fraction of TiCl4 consumed after rapid compression in the RCM, plotted against the peak temperature reached. The lower limit shows a simulation with tc ¼10 ms, the upper limit 40 ms, and the filled diamond the tc that was reported for the corresponding measurement (white circle) from [10].
1768
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770
our simulations than the reported measurements, but, although no error bars were reported for the measurements, we expect the uncertainties in the temperature estimations to be significant (see Section 5).
higher than that observed in [5]. This may be partially because mass transport limitations were not modelled. Also, Ghoshtaghore’s rate expression is for much higher O2 partial pressures, so would be expected to over-predict wall reactions in the conditions modelled here.
4.2. Plug Flow Reactor (PFR) 4.3. Particle properties The combined simulation coupling the gas-phase kinetic model to the particle dynamics model was used to simulate the tubular flow reactor used by Pratsinis et al. [5] to measure the overall reaction kinetics. In [5], a flow of 0.2 mol% TiCl4 and 1 mol% O2 in Argon was fed to a 3.18 mm-I.D. tube, a 30 cm section of which was heated to 973–1273 K in a furnace. The concentration of TiCl4 leaving the flow tube (C o ) was measured using FTIR spectroscopy. Assuming that the reaction rate is first order in [TiCl4] with Arrhenius kinetics (keff ¼ A expðEa =RTÞ) and that all reaction occurs during the time t that the gas spends in an isothermal zone at temperature T then
ðln ðC o =C i ÞÞ=t ¼ keff ¼ A expðEa =RTÞ:
ð3Þ
A plot of ln ððlnðC o =C i ÞÞ=t Þ versus 1=T (Fig. 5) would thus give a straight line from which the apparent activation energy Ea and pre-exponential constant A can be determined. This is the analysis performed in [5], and our simulation results are processed the same way for easy comparison. The assumption of constant T for a time t was only made for the analysis: the simulations covered the entire temperature profile (given in [5]) and included thermal expansion of the gas as the temperature varied. The simulations using our new gas-phase kinetic model (Fig. 5, solid lines) consume much less TiCl4 than the observations reported in [5]. At 1133 K the keff fitted to the simulation is between 14 and 45 times lower (slower rate) than that fitted to the experiment, depending on which residence time is simulated. The apparent activation energy is higher in the simulation, so at 993 K the expt sim discrepancy increases to keff =keff ¼ 5:2 103 in the worst case. Possible reasons for this discrepancy will be discussed in a later section. The simulations were repeated allowing for deposition and reaction of TiCl4 on the inside surface of the tube, using the expression for surface growth of a TiO2 film given by Ghoshtagore [15] (Eq. (2)). In these simulations (Fig. 5, dotted lines), 500 times as much TiCl4 was consumed through reaction on the walls as through reactions in the gas phase, and the overall consumption of TiCl4 was
The flow tube reactor was simulated using the coupled population balance model. This allows access to a wealth of information far beyond the rate of consumption of TiCl4. Fig. 6 shows the size distribution of titania particles at the end of the simulation (40 cm downstream from the midpoint of the tube, when the temperature has returned to 298 K) and Fig. 7 shows TEM-style images of some of the particles. These simulated TEMs are created from the list of primary particles that is produced by the population balance code as per the method described in [36]. As the temperature of the PFR is increased, so does the average size of the particles produced. This is the same as the trend observed by [37] in a similar PFR. 4.4. Comparison of rates Although the complex pressure-dependent chemistry of this system cannot be fully described by a first-order rate expression it can still be interesting to compare the different measurements and predictions in terms of an apparent first-order rate coefficient. The determination of an effective rate coefficient from the plug flow reactor experiments was described in Section 4.2. The best fit Arrhenius expressions reported in [5] for four different O2/TiCl4 ratios are plotted on an Arrhenius plot in Fig. 8. Comparable rate coefficients for the rapid compression machine and the detailed kinetic model were found as follows. 4.4.1. From the detailed kinetic model A batch reactor at constant temperature and pressure and initial composition of 1 mol% TiCl4 and 5 mol% O2 in Argon was simulated for a time t ¼ minðthalf ; 1sÞ, which is the lesser of the time at which C o =C i ¼ 0:5 or 1 s. The choice of maximum reaction time (1 s) has only a small effect on the results. The apparent first-order rate coefficient keff was determined according to
keff ¼ ðln ðC o =C i ÞÞ=t:
ð4Þ
The simulations were repeated at pressures from 0.1 to 100 bar to investigate the pressure-dependence of the model. As shown in
4 2
100 T=993 K T=1073 K 10
−2 d P/ d log Dp
ln(k eff )
0
−4 −6
0.49 s
−8
0.75 s
T=1133 K T=1293 K
1
0.1
1.1 s −10 0.85
0.9
0.95
1.0
103 K/T Fig. 5. Arrhenius plot of plug flow reactor at three different residence times. Unconnected points are measurements reported in [5], lines are our simulations: dotted with deposition on the reactor walls and solid without.
0.01 0.0001
0.001
0.01 Dp (µm)
0.1
1
Fig. 6. Particle size distributions from the PFR simulation with residence time in the hot-zone of 0.49 s.
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770
1769
Fig. 7. Simulated TEM Images from the exit of the PFR simulation with residence time in the hot-zone of 0.49 s and furnace temperatures 993, 1017, 1073 and 1293 K (left to right).
Fig. 8, the apparent rate constant is much higher at higher pressures, increasing almost linearly with pressure. 4.4.2. From Raghavan’s RCM data To determine an effective first-order rate coefficient from the RCM measurements requires some assumptions. First, that the temperature follows a Lorentzian with HMHM of t c and peak temperature of T p . Second, that the temperature-dependence of the rate coefficient follows expðEa =RTÞ. Ea was taken to be 108 kJ/ mol as given by Raghavan, but the results are not very sensitive to this choice. With these assumptions the effective rate coefficient was determined from a single compression experiment according to
keff ¼
! ln ððN=N0 f Þ=ð1 f ÞÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : tc pRT p =Ea
ð5Þ
Where f is the fraction of the gas that is compressed isothermally according to the hot-core model. There are some difficulties with determining keff from the observed ðN=N 0 Þ values: first, if N=N0 P 1 then the apparent rate is negative; second, if N=N 0 6 f then the apparent rate is greater than infinity. The uncertainties of keff values from experiments with high and low conversions are therefore very large, making lnðkeff Þ positive and negative infinity, respectively. To illustrate this the bars in Fig. 8 shows the error in lnðkeff Þ that would result from an error of 0.1 in the measurement of ðN=N0 Þ. 5. Discussion A rapid compression machine seems to be a useful tool for measuring gas-phase kinetics at high temperatures and pressures, but at less expense than a shock-tube. However, uncertainties in the temperature measurements are likely to be significant. The peak temperature T p was deduced from measurements of pressure, volume, and the starting temperature T 0 , giving three causes for con-
Our Simulation 10
Raghavan RCM Pratsinis PFR
ln(k)
5 0 -5
100 bar 10 bar 1 bar 0.1 bar
-10 0.5
1.0
1.5
cern. First, an error in measuring T0 is amplified in T p . Second, the deduced T p is very sensitive to errors in measuring the maximum piston displacement. Third, the assumptions necessary for this deduction are hard to verify. Unfortunately, Raghavan does not give estimates of the errors in N=N 0 or T p in Fig. 4, but we expect they are not small enough to justify adjusting parameters in our kinetic model to fit these experiments. The 20% dead-zone in the model ensures that no simulation will ever match the observations in which ð1 N=N 0 Þ > 0:8, whichever kinetics are used. Plug flow reactors are often used to measure rates, and the careful measurements in [5] seem like a useful check for any kinetic model of TiCl4 oxidation. At first glance the disagreement of four orders of magnitude between the first-order rate fitted to the experiments and that fitted to the simulations suggests a problem with the new gas-phase kinetic model used in the simulations. Indeed, this may be the case, and further work on the model is recommended, especially the initiating reactions. However, there is an alternative possibility. Due to the small internal diameter of the tube, the surface-area to volume ratio is high, and including surface reaction on the tube wall has a large effect on the simulations. It is possible that some of the observed consumption of TiCl4 occurred on the walls of the tube, rather than in the gas phase. The rate expression derived from these measurements, however, has been widely used to simulate gas-phase reactions, and has often been used as an inception rate for TiO2 particle population balance models. Yoon et al. [38] observed that, in their tubular reactor of diameter 2.7 cm, up to 92% of TiO2 production was on the walls; although a different experiment, this illustrates that in at least some conditions, wall reactions can dominate. Further investigation of both gas and surface kinetics may be required to resolve this issue. Experimental investigations into the influence of tube diameter on the overall rate of decomposition could be particularly instructive. The uncertainty in the kinetic model is difficult to assess fully and a proper treatment is beyond the scope of this paper. Spectral expansion techniques such as those in Sheen et al. [39], Russi et al. [40] and Xiu et al. [41] can be used to rigorously quantify uncertainties in the kinetic rate parameters. Further work using techniques like this will enable the determination of the overall errors associated with the model. As shown in Fig. 8 the detailed kinetic model predicts that the overall reaction rate is pressure-dependent, and the temperaturedependence is too complicated to describe with a single Arrhenius expression. This decreases the accuracy of rate expressions determined experimentally over a narrow range of temperature and pressure, especially when extrapolated beyond that range, for example to model industrial reactor conditions. The current detailed model is much more versatile in terms of operating conditions. 6. Conclusions
(1000 K) / T Fig. 8. Arrhenius plots of our model at different pressures with Pratsinis [5] and Raghavan [10] for comparison.
The first detailed gas-phase kinetic model for TiCl4 oxidation [11] has been improved using DFT and VTST calculations, which
1770
R.H. West et al. / Combustion and Flame 156 (2009) 1764–1770
have provided new thermochemical and kinetic data. This improved kinetic model is in agreement with the experimental data from a rapid compression machine [10], with the results of the model falling within the estimated uncertainties of the experimental measurements. Simulations of a plug flow reactor, however, suggest that surface reactions on the reactor walls may have some impact on the experimental measurements [5]. Many models simulating TiO2 formation from the gas phase have relied on rate expressions derived from these experimental measurements. The improved kinetic model thus raises an important question in the field of modelling TiO2 particle formation. Finally, the improved gas-phase kinetic model was coupled to a particle population balance model (PBM) to investigate the influence of PFR temperature on particle size distributions. The traditional surface-volume PBM was extended to track the growth and sintering of primary particles within each particle. This allowed creation of TEM-style images of the agglomerate particles from the simulation, which are presented herein. Acknowledgments The authors would like to thank Prof. Stephen Klippenstein for help with the Variflex calculations, Tioxide Europe Limited (TEL) for the financial support of R.A.S., the EPSRC for Grant No. EP/ E01724X/1, and both the EPSRC and TEL for the financial support of R.H.W. through an Industrial CASE studentship. Appendix A. Supplementary data
[13] [14] [15] [16] [17]
[18]
[19] [20] [21] [22] [23] [24]
[25]
[26]
Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.combustflame.2009.04.011. [27]
References [1] J. Emsley, Molecules at an Exhibition, Oxford University Press, 1999. [2] R.A. Gonzalez, C.D. Musick, J.N. Tilton, Process for controlling agglomeration in the manufacture of TiO2, US Patent 5,508,015, 1996. [3] J.C. Deberry, M. Robinson, M. Pomponi, A.J. Beach, Y. Xiong, K. Akhtar, Controlled vapor phase oxidation of titanium tetrachloride to manufacture titanium dioxide, US Patent 6,387,347, 2002. [4] R.H. West, G.J.O. Beran, W.H. Green, M. Kraft, J. Phys. Chem. A 111 (18) (2007) 3560–3565. [5] S.E. Pratsinis, H. Bai, P. Biswas, M. Frenklach, S.V.R. Mastrangelo, J. Am. Ceram. Soc. 73 (1990) 2158–2162. [6] P. Koukkari, J. Niemela, Comput. Chem. Eng. 21 (1997) 245–253. [7] P.T. Spicer, O. Chaoul, S. Tsantilis, S.E. Pratsinis, J. Aerosol Sci. 33 (2002) 17–34. [8] S. Tsantilis, S.E. Pratsinis, J. Aerosol. Sci. 35 (2004) 405–420. [9] N.M. Morgan, C.G. Wells, M.J. Goodson, M. Kraft, W. Wagner, J. Comput. Phys. 211 (2) (2006) 638–658. [10] R. Raghavan, Measurement of the High-Temperature Kinetics of Titanium Tetrachloride (TiCl4) Reactions in a Rapid Compression Machine, Ph.D. thesis, Case Western Reserve University, 2001. Available from:
[28] [29] [30]
[31] [32] [33] [34] [35] [36]
[37] [38] [39] [40] [41]
O. Inderwildi, M. Kraft, ChemPhysChem 8 (2007) 444–451. S.E. Pratsinis, P.T. Spicer, Chem. Eng. Sci. 53 (1998) 1861–1868. R.N. Ghoshtagore, J. Electrochem. Soc. 117 (1970) 529–534. J.A. Montgomery Jr., M.J. Frisch, J.W. Ochterski, G.A. Petersson, J. Chem. Phys. 112 (15) (2000) 6532–6542. M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, J.A. Montgomery, Jr., T. Vreven, K.N. Kudin, J.C. Burant, J.M. Millam, S.S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G.A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J.E. Knox, H.P. Hratchian, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, P.Y. Ayala, K. Morokuma, G.A. Voth, P. Salvador, J.J. Dannenberg, V.G. Zakrzewski, S. Dapprich, A.D. Daniels, M.C. Strain, O. Farkas, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J.V. Ortiz, Q. Cui, A.G. Baboul, S. Clifford, J. Cioslowski, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, C. Gonzalez, J.A. Pople, Gaussian 03, Revision C.02, Gaussian, Inc., Wallingford, CT, 2004 (2003). S.J. Klippenstein, A.F. Wagner, R.C. Dunbar, D.M. Wardlaw, S.H. Robertson, J.A. Miller, VariFlex, software package, 2002. Available from: