J. Phys. Chem. 1996, 100, 9323-9333
9323
Detailed Numerical Simulations in Flow Reactors: A New Approach in Measuring Absolute Rate Constants J. Segatz McKinsey & Company, Taunusanlage 21, 60325 Frankfurt, Germany
R. Rannacher Institut fu¨ r Angewandte Mathematik, UniVersita¨ t Heidelberg, INF 293, 69120 Heidelberg, Germany
J. Wichmann, C. Orlemann, T. Dreier,* and J. Wolfrum Physikalisch Chemisches Institut, UniVersita¨ t Heidelberg, INF 253, 69120 Heidelberg, Germany ReceiVed: December 4, 1995; In Final Form: February 23, 1996X
A detailed investigation is presented to simulate the reactive flow field in a low-pressure flow reactor for kinetic studies. This has been done to improve existing methods for evaluating data from isothermal flow kinetic measurements. The full Navier-Stokes equations for compressible flows including transport phenomena and chemical reactions have been solved numerically for the low Mach number case. By using splitting techniques for variables and spatial dimensions, the calculation time could be reduced by more than 2 orders of magnitude. This reduction allows the repeated application of the solver to adjust parameters in the kinetic model organized as an optimization problem and give best agreement between experiment and calculation. The model results for the nonreactive flow field have been verified by comparison to imaging measurements via two-dimensional laser-induced fluorescence of acetone tracer gases for the visualization of diffusive mixing. Numerical results of full reactive flow simulation have been compared with the measurement of elementary relaxation processes and vibrational energy transfer in collisions of vibrationally excited hydrogen and deuterium molecules. Spatially resolved axial and radial concentration profiles of both species were obtained at room temperature using coherent anti-Stokes Raman spectroscopy (CARS). From the detailed numerical simulation evaluated wall deactivation probabilities at 300 K for H2(V)1) f(wall) H2(V)0) (Ia) of γw ) (1.5 ( 0.3) × 10-3 s and thermal rate constants for vibrational energy transfer H2(V)1) + D2(V)0) f H2(V)0) + D2(V)1) (IIa) of kvv ) (6 ( 0.5) × 109 cm3 mol-1 s-1 were derived using an optimization procedure specially adapted to the present kinetic problem. They are larger, respectively, by a factor of 2 (wall deactivation) and 1.4 (vibrational energy transfer) compared with values obtained from a plug-flow evaluation. While the kvv data from different experiments are now in excellent agreement, theoretical results using recent ab initio potentials still differ by a factor of 2.
Introduction Since the beginning of kinetic measurements in dischargeflow systems in the 1920s following the pioneering work of Wood1 and Bonhoeffer2 the flow tube technique kept its importance in modern experiments as one of the most powerful tools for the determination of elementary chemical thermal reaction rate constants. Whereas new technologies, especially in spectroscopic detection methods, have extended its range of applicability, the basic principle of flow tubes remains the same: mixing of reactants takes place upstream in a mixing section and consumption of reactants or buildup of products is followed along a measurement section of the flow tube by some detection method for atoms, radicals, or molecules to deduce a rate constant from measured axial concentration profiles. An assumed mean flow velocity converts the axial coordinate (distance between the point of mixing and detection) into reaction time. The reaction rate constants of interest can then be deduced by modeling the homogeneous chemical reaction system. For most carefully designed kinetic measurements this “plug flow” assumption is a useful approximation and has produced a wealth of trustworthy kinetic data. However, the method is known to bear systematic errors, since it is based on X
Abstract published in AdVance ACS Abstracts, April 15, 1996.
S0022-3654(95)03578-7 CCC: $12.00
the approximation of a perfect decoupling of chemical and hydrodynamic processes in the flow tube. Especially in the mixing section of the reactor this assumption is not valid. Concentration profiles of reactants and products in a flow tube usually show spatial variations in radial direction caused by the characteristics of the underlying flow field. Generally, the flow velocity is higher in the central region of the reactor than near the wall, giving rise to different residence times of the reactant species depending on their radial position. Those travelling in low speed domains near the wall experience a longer reaction time than those travelling near the center at the same axial position. In addition, for unstable species such as radicals or internally excited atoms and molecules, despite diffusional transport heterogeneous reactions and relaxation on reactive surfaces also entail radial concentration gradients depending on geometric constraints of the flow tube and wall activity. Often concentrations are deduced from spatially integrating optical detection methods, i.e., absorption, which disregard these spatial variations. In order to favor diffusive processes, which minimize radial concentration gradients, a flow tube traditionally is operated at low pressure (less than 1.3 kPa), in the so-called “plug flow limit”.3 On the other hand a rapid radial diffusional transport will increase the relative importance of wall reactions with respect to the gas phase reactions of © 1996 American Chemical Society
9324 J. Phys. Chem., Vol. 100, No. 22, 1996
Segatz et al.
interest. Besides viscous pressure drop, back diffusion will become important: a diffusional transport process in the axial direction governed by axial concentration gradients that are caused by the chemical reactions. Limitations based on detector sensitivity impose further restrictions upon operating at low pressures. Therefore, for the special case studied in this work, where relatively slow flow velocities and the constraint for measuring in a very early stage of mixing exist, it is difficult to achieve true plug flow conditions. The present study is intended to permit numerical modeling studies to analyze such systems in greater detail than was previously available. In the past, for the evaluation of kinetic rate constants many studies have been conducted to consider the various mentioned aspects by approximations to the ideal reactive flow. One of the first detailed discussions was given by Kaufmann.4 He addressed the basic problems arising from viscous flow and diffusion and derived approximate expressions for axial diffusion as well as for the treatment of wall reactions in terms of a pseudo homogeneous rate constant and an estimation of the radial concentration gradient due to homogeneous and heterogeneous reactions in a Poiseuille flow. There have been various attempts to extend the flow tube technique to pressures higher than imposed by the plug flow limit, by solving the two-dimensional continuity equation for a single reactant species in chemical first- or second-order reaction:
(
)
∂2n 1 ∂ 1 ∂n ∂n + kn )D 2+ ∂z r ∂r r ∂r ∂z
Vz(r)
(1)
with Vz(r) being the flow velocity at axial and radial position z and r, respectively, D the diffusion coefficient, n the reactant concentration, and k the first-order rate constant. Poirier and Carr5 were able to correct observed axial concentration decays up to total pressures of 2 kPa helium by numerically solving eq 1 for a first- or second-order homogeneous reaction and a first-order heterogeneous reaction. Their assumptions were a fully developed laminar flow, a top hat initial concentration profile, and negligible axial diffusion using experimentally determined wall loss and diffusion constants. Ogren6 developed an approximate analysis that considers homogeneous and heterogeneous first-order reactions neglecting axial diffusion, and assuming fully developed laminar flow. Westbrook et al.7 developed a technique for studying kinetics under turbulent, rather than conventional laminar flow conditions. They report results for studies of high-temperature combustion processes at atmospheric pressures. Keyser8 implemented an analytic solution of the continuity equation, initially found by Walker9 and later developed by Brown,10 into the data evaluation procedure, again assuming Poiseuille flow. Their results were valid up to 10 kPa of helium. Abbat et al.11 used experimentally determined velocity profiles as well as measured axial/radial reactant concentration profiles to solve eq 1 numerically under pseudo-first-order conditions without the assumption of fully developed laminar flow. They were able to extend the range of applicability up to 50.8 kPa. Seeley et al.12 applied a turbulent flow technique to the evaluation of elementary reaction rate constants. For Reynold numbers larger than 2000, the flow separates into a turbulent core with constant velocity and a laminar sublayer near the wall moving much slower. The authors showed that, as long as detection is limited to the core, the flow can be well described by a one-dimensional continuity equation using the measured core velocity and an effective diffusion coefficient, which incorporates mixing by both molecular diffusion and turbulent processes. Nonetheless, so far a reliable and complete treatment of the mixing region has not been done, so that important information
from the initial stage of reaction remains unconsidered. Also, within the previous approaches it is often difficult to implement a complex kinetic mechanism. This is especially important, since oftentimes measurements have to be carried out under experimental conditions which do not allow a description in terms of a one species continuity equation such as (1). Therefore, in a more reliable evaluation of rate constants from experimental data, it is desirable to take into account all relevant physical and chemical processes occuring in reactive flows, without having to depend on simplifying mathematical assumptions. This can only be done by a detailed numerical treatment. A practicable simulation again places demands on the experimental layout, especially the design of the flow reactor. The purpose of the present work is to show that the detailed modeling of reactive flow fields within a suitable reactor for kinetic studies is a powerful tool for the experimental determination of elementary rate constants. Comprehensive analytical and experimental validations guarantee that the simulations properly describe transport and chemical processes in the reactive flow. It will be shown that, even for low pressures, the approximations made within the simplified plug flow evaluation can lead to underestimation of rate constants. As a model system for a simple kinetic process the heterogeneous relaxation of vibrationally excited hydrogen (H2(V′′)1)) and its energy transfer in collisions with deuterium (D2(V′′)0)) is considered: wall
H2(V′′)1) 98 H2(V′′)0) H2(V′′)1) + D2(V′′)0) f H2(V′′)0) + D2(V′′)1)
(Ia) (IIa)
Reaction IIa is of interest because of its model-like character as a simple example for the theoretical understanding of vibrational energy transfer (VV) processes in molecular collisions. There is still a discrepancy between results of theoretical approaches and experimentally determined rate constants.13 Reactions Ia and IIa also are of importance in the context of experimental rate measurements for the vibrationally mediated hydrogen exchange reaction and its isotopic variants, the prototype of a bimolecular reaction in the gas phase.14 In this reaction the presence of vibrationally excited reactants requires the knowledge of all deactivation pathways, such as heterogeneous relaxation and vibrational energy transfer, which compete with the reactive channel. Experimental Section Flow System. Gas handling and the principal outline of the flow system used in the present investigations are very similar to the setup described in ref 15. Briefly, the reactant gases (H2, purity 99,999%, D2 (99,7%)) and buffer gas (He, 99,999%) are metered by flow controllers (Tylan) and further purified in liquid nitrogen cold traps before entering the flow reactor. It consists of the reactor head where gases enter the actual flow tube. Both parts can be attached to one another by an O-ring sealed flange system which enables an easy change of all design parameters of the flow reactor. We applied two kinds of reactor heads. The conventional one consisted of a T-shaped Pyrex glass tubing (Figure 1a) where the gas flows impinge on each other directly and cause a fast, turbulent mixing. In the literature, in most kinetic studies using flow tubes measurements are made further downstream from the inlet ports where mixing is complete, and close to plug flow conditions are established. Due to the complex flow conditions, concentration measurements in this mixing section cannot be used for our data evaluation. Therefore, a second reactor head was designed (Figure 1b) in order to circumvent these difficulties in the data evaluation
Numerical Simulations in Flow Reactors
Figure 1. Design of the Pyrex/quartz mixing heads attached on top of the flow tube: (a) turbulent mixing, (b) generation of a laminar flow through a fourfold gas inlet symmetrically arranged around the central tube. For details see text. The geometry is as follows: Le ) 45 mm, do ) 32 mm, di ) 10 mm, da ) 14 mm, dB ) 2 mm, aB ) 12 mm (5×), 15 mm (5×), 20 mm (6×).
and simultaneously provide a swirl-free, two-dimensional rotationally symmetric flow field, which is accessible to a detailed numerical treatment from the beginning of the reaction. One reactant is admitted through an inner axial tube (diameter 10 mm, wall thickness 2 mm), the other through four equally spaced side arms. Reactants get in contact at the outlet of the central inlet tubing, the length of which is large enough to guarantee fully developed laminar flow fields in both the inner and outer gas flow. Mixing is mainly governed by diffusive transport. For comparison all kinetic measurements have been performed under the same conditions in both the laminar and the conventional flow system, and have been evaluated either by numerical simulations or under the plug flow assumption, respectively. The main flow tube (see Figure 2) consists of a straight, 32 mm internal diameter section equipped with an array of diametrically opposed 2 mm diameter holes in the wall to allow for optical CARS diagnostics with focused laser beams. In order to avoid damage of window material located too close to the beam foci, small diameter 70 mm long tubes fitted with quartz windows at one end are glued to the central flow tube so that beam intensities are low enough on the window surfaces. This precaution also reduces the generation of spurious nonresonant CARS signals within the window material. Various materials and coatings of the flow tube wall have been tested to minimize heterogeneous relaxation of vibrationally excited species. From CARS measurements it was found that stainless steel coated with Teflon (Du Pont FEP 856200) leads to the lowest wall deactivation probability for H2(V′′)1), 2 times lower than the value obtained with HF-cleaned Pyrex glass.15 As is shown in Figure 2, the flow reactor was positioned into the beam path and aligned at each axial measurement location by a translation/ rotation stage. Low pressure is maintained by a mechanical pump (Alcatel, 2063), and pressure can be read from a precision capacitance manometer (MKS Baratron, 290H). Vibrationally excited hydrogen molecules were generated in two microwave cavities assembled to the side arms of the respective reactor
J. Phys. Chem., Vol. 100, No. 22, 1996 9325
Figure 2. Three-dimensional view of the flow reactor used for kinetic measurements. For the generation of vibrationally excited species microwave discharge cavities are attached to the quartz inlet tubing ((2) in figure 1). Shown are the window extension tubes for transmitting the laser beams. The reactor is attached to translation and rotation stages to allow for proper spatial alignment.
Figure 3. Schematic of two-dimensional laser-induced fluorescence experimental setup for the monitoring of the spatial distribution of acetone tracer species inside the flow tube. As depicted, the light sheet excites a horizontal cross section of the flow downstream the inlet manifold and detection is through a quartz window at the bottom next to the pumping port.
head (cf. Figure 2), and driven by a magnetron microwave generator (2450 GHz) at powers between 20 and 100 W. 2D LIF Setup. Rotational symmetry and the absence of vortices and swirls in the flow were checked experimentally by adding small amounts of acetone as a tracer to the outer or inner gas flow and visualizing their spatial distribution by planar laser induced fluorescence (2D-LIF) techniques16 taken in vertical and horizontal sections at different axial locations along the flow tube (see Figure 3). For excitation the fourth harmonic of a pulsed Nd:YAG laser at 266 nm is formed into a light
9326 J. Phys. Chem., Vol. 100, No. 22, 1996
Figure 4. CARS experimental setup: λ/4 and λ/2: quarter- and halfwave plates; SEV, photomultiplier tubes.
sheet by a cylindrical lens (f ) 700 mm). Care was taken to place the focal line behind the flow reactor which ensured a (slowly converging) sheet thickness of 500 µm in the measurement section. For these experiments the regular flow tube was replaced by a Suprasil tube of equal dimensions. Planar LIF concentration maps of the tracer gas were imaged perpendicular to the plane of excitation either through the side walls (vertical cuts) or, as is shown in Figure 3, through a quartz window at the bottom of the flow tube (horizontal cuts). The broadband fluorescence signal was detected through a blue glass filter between 350 and 550 nm by a CCD camera (Proxitronic Nanoscan PC 1811) with adjustable image intensifier. For a detailed description of 2D-LIF experiments for concentration measurements we refer to refs 16 and 17. CARS setup. The CARS experimental setup is depicted in Figure 4. A pulsed (8 ns) single mode Nd:YAG-laser (Spectra Physics, GCR-3) and a narrowband (∆ω ) 0.1 cm-1) dye laser (LUMONIX, HD500) pumped by the same laser provide pump and Stokes beams for generating CARS signals in hydrogen and other gases in a collinear beam geometry similar to a setup described in ref 15. Pyridin I dye (660-740 nm) and DCM (610-675 nm) in methanol provide the appropriate Stokes wavelengths for hydrogen and deuterium, respectively. Galilei telescopes in each beam path enable the separate adjustment of the beam waists for proper focusing and recollimation of both beams to a common focal spot using f ) 300 mm focal length lenses. For the measurement of radial profiles inside the reactor these can be translated along the optical axis to move the focal spot along the beam direction. The CARS signal emerging from the focal volume is separated from pump and Stokes beams by three longwave pass dichroic mirrors, interference filters and a 0.25 m monochromator (McPherson, 218) before it is detected by a photomultiplier (EMI 9635B). The captured and integrated (SR250) electrical signals are plotted on a pen recorder. To eliminate systematic errors in concentration measurements due to interference of nonresonant signals from buffer gases and windows with the resonant part of the third-order susceptibility of the detected species, the nonresonant contribution to the CARS signal was rejected by polarization selective excitation and detection.18 This technique reduced the nonresonant background by 3 orders of magnitude, giving at 5.33 mbar (5.3 × 102 Pa) total pressure a detectivity of 0.06% corresponding to a quantum state specific detection limit of 7 × 1012 molecules/ cm3. Measurement Procedure and Data Evaluation. Vibrationally excited hydrogen molecules H2(V′′)1) were generated
Segatz et al. by microwave discharges in the sidearms of the mixing head (cf. Figure 2). Populations in higher vibrational states were below the current CARS detection limit. Dilution ratios in rare gas helium and discharge conditions are optimized for minimum atom yield and maximum vibrational excitation. When the microwave discharge was switched on, within the measurement error, there was no significant (e.g., by CARS) temperature variation detectable down the flow tube axis. Hydrogen atom concentrations were determined by quantitative gas phase titration with ethylene using small quantities of NO to generate a detectable chemiluminescence from HNO*.15,19,20 Generally, for highest yield in [H2(V′′)1)] with simultaneous low atom yield in the discharge a high (almost undiluted) hydrogen partial pressure and a low microwave power (20 W) in the two side arms were preferred. Typical atom partial pressures were in the range of 2 Pa, which corresponds to a dissociation yield of 0.5%. Hydrogen molecules, on the other hand, can be excited vibrationally with an efficiency of 0.8%, which gives a concentration of 3.5 × 1014 cm-3 at the present conditions. The CARS lines of individual Q-branch transitions of hydrogen are well separated and their full width at halfmaximum (fwhm) at room temperature and low pressure are entirely determined by Doppler broadening. Due to efficient suppression of the nonresonant background, the measured CARS signal intensity is proportional to the square of the population difference of the Raman coupled states. For determination of relative population densities from the measured population differences the assumption is made that for the highest vibrational level probed with a detectable signal intensity the population in the upper Raman level is zero. Thus, the relative population densities of all lower levels can be calculated successively, if account is taken of the different Raman cross sections of the respective transitions. Conversion to absolute concentrations was possible from the known volumetric flow rates and total pressure (533 Pa). Generally, the most intense Q1 and Q2 CARS lines for hydrogen and deuterium, respectively, are used for concentration measurements. When vibrational population ratios are determined, care was taken to employ low laser beam energies (Ep < 1 mJ, ES < 0.2 mJ) to exclude errors from vibrational quantum number dependent saturation.21 This was verified by monitoring the CARS signal intensity as a function of pump and Stokes beam energy over several orders of magnitude (for beam energies between 2 µJ and 6 mJ). Only below the above-cited energy levels was there a linear and quadratic increase of the signal intensity from the Stokes and pump beam energy, respectively. Modeling and Simulation of the Transport-Reaction Problem In a flow reactor several interacting processes like convection, diffusion, gas phase and wall reactions take place simultaneously. Due to this interaction experimental data are often difficult to interpret. By means of a detailed numerical simulation, these various aspects can be distinguished and the interacting processes of mixing and reaction are easier to understand.22 Governing Equations. The intention for a numerical simulation is to provide two-dimensional profiles for concentration, temperature, density and velocity fields. Therefore, the governing equations are of Navier-Stokes-type combined with species and mass conservation:23
mass:
∂F/∂t ) -∇‚(Fv)
(2)
species:
∂Fi/∂t ) -(∇‚Fiv) - (∇‚jd,i) + ω˘ iMi
(3)
Numerical Simulations in Flow Reactors
J. Phys. Chem., Vol. 100, No. 22, 1996 9327
momentum: ∂Fv/∂t ) -[∇‚Fvv] + [∇‚(-pI + τ)] + Fg (4) energy:
∂FU/∂t ) -(∇‚FUW) - (∇‚jq) - p(∇‚v) + (τX∇v) + q˘ (5)
with F, Fi, and v the total density, species density, and velocity vector, respectively. Further, jd,i, ωi, Mi, and p are the diffusion flux, the chemical source term, specific mole mass, and pressure, respectively. The stress tensor, gravity vector, internal energy, heat flux, and energy source terms are denoted by τ, g, U, jq, and q, respectively. Furthermore, the following relations hold: 3,23
R p)F T M h
ideal gas law: stress tensor:
(6)
2 τ ) - µ(∇‚v)I + µ[∇V + (∇v)T] 3
(7)
wi jd,i ) jDd,i + jTd,i ) -F DDi∇xi xi
diffusion flux:
DTi ∇T (8) T ns
heat flux:
jq ) j q + j q ) -λ∇T + ∑hijd,i c
d
(9)
i
detailed chemistry:
nr
ns
l)1
i
ω˘ ) ∑{(a˜ il - ail)kl∏ciail}
{ }
kl ) AlTβl exp
-Eal RT
(10)
(11)
with R universal gas constant, M mean molar mass, T temperature, µ viscosity, DD mass diffusion coefficient, DT thermal diffusion coefficient, wi mass fraction, xi mole fraction, λ thermal conductivity, hi specific enthalpy, ail stoichiometric number, kl reaction rate constant, ci species concentration, Al preexponential factor, β1 temperature exponent, Eal activation energy, ns number of species, and nr number of reactions. Detailed diffusion and heat transport as well as specific gas phase and surface reaction mechanisms are taken into account. Vibrationally excited species are treated “state-selective”, i.e., each species with a different vibrational level represents a new species with specific thermodynamic and physical properties, for which enthalpies and transport coefficients are changing correspondingly. Initially, close to 200 reaction steps were considered for the present system. However, a sensitivity analysis using a zerodimensional kinetic algorithm revealed that many of them were of marginal importance. Therefore, to handle the present full reactive flow problem with a reasonable computational speed, including the optimization procedure, the reaction mechanism was reduced to the sequence listed in Table 1. Because of the flow tending to the incompressible case (Mach number < 0.01), special efforts are made to stabilize the variable-density algorithm also for this kind of flow.32 Here the implementation of a pressure-correction equation33 proved to be most efficient, since it eliminates the acoustic wave propagation and therefore reduces the numerical stiffness of the system. Instead of using the nonstationary continuity equation for evaluating the density distribution (taking the pressure from the ideal gas law, which has a very weak coupling in the incompressible limit), it is more reasonable to obtain the pressure
field from a pressure correction equation of elliptic type (which is generated by combining continuity and momentum equations) and calculating the density by means of the ideal gas law.34 The resulting equation for the pressure field reflects the nearly elliptic behavior of the propagation of pressure disturbances, since the Mach number of an incompressible flow tends to be infinity and so the pressure-velocity coupling is given everywhere in the domain. Spatial Discretization. For spatial discretization a finitedifference scheme on an adaptive tensor product grid is used. Regions of steep gradients in the computational domain are resolved with higher resolution controlled by an error estimate. In this case a refinement strategy is used, which allows a limitation of the discretization error. In a first step a stationary solution on a coarse grid is computed (level 0 solution). The grid then is globally refined by taking half the grid spacings, and again the steady state is calculated (level 1 solution). By comparing level 0 and 1, zones of changes in the solution can be identified which are selected for further local refinement. As a criterion, whether an element has to be refined or not, the relative change of the solution is compared with a given tolerance limit. Proceeding in this manner, the grid is locally refined, until the “response” of the solution (e.g., the relative change) is smaller than the given tolerance limit everywhere in the computational domain.35 Differential operators are approximated with a centered ninepoint stencil. The geometry is represented in two-dimensional (axial and radial) cylindrical coordinates. A “staggered-grid” technique was used to obtain oscillation-free pressure and velocity fields, where different control points are used for scalar variables and velocities. Temperature, density, and species composition are calculated at central grid points, whereas the momentum equations are considered at points that are shifted half a grid spacing in x and y direction for the x and y velocity components, respectively. For stabilizing convective terms, upwind discretization of first order is optionally available in the code.32 Boundary Conditions. For reasonable numerical results the implementation of correct boundary conditions is important.34 Four different types of conditions are applied in this simulation: inflow conditions are of Dirichlet (fixed values) type for all variables; only density is calculated by a time-differential equation. For investigating the premixing effect of backdiffusion on the inflow concentration, the boundary of the computational domain was shifted upstream. The flow tube axis is characterized by the symmetry condition of vanishing radial gradients for all variables. Only the radial velocity is set to zero. The tube wall is characterized by no-slip conditions for all velocities (i.e., velocities are set to zero). Density is calculated by time-differential equations, and equilibrium is assumed of the diffusive flux and source term for temperature and species concentrations. The outflow boundary is treated by one-sided differential equation for all variables (for axial direction only first-order derivatives are calculated). Time Integration. Although the flow problem is of stationary type, solving the above equations by nonstationary or quasistationary methods is of advantage. The best efficiency is obtained by using a splitting scheme for variables and spatial dimensions. Splitting of the spatial directions is performed by approximate factorization (AF), an ADI (alternating direction implicit) similar method. Especially for flows showing a dominant flow direction without reverse-flow areas, this method is well suited. Such conditions are fairly fulfilled in the rotationally symmetric geometry of the present situation. When using the AF method, the solution of the system of equations splits up into two steps, where each forces the coupling into
9328 J. Phys. Chem., Vol. 100, No. 22, 1996
Segatz et al.
TABLE 1: Thermal Rate Constants k ([m3/(mol s)]) and Deactivation Probabilities γ Used in the Kinetic Modeling for Vibrational Energy Exchange and Relaxation of Hydrogen Isotopes in the Flow Reactor no. Ia Ib Ic
reaction H2(V′′)1) D2(V′′)1) H+H
Heterogeneous Relaxations wall H2(V′′)0) 98 wall D2(V′′)0) 98 wall H2(V′′)0) 98
γ(300 K)
ref
1.49(-03)a 1.49(-03) 1.00(-04)
this work this work 31
k(300 K)
a
IIa IIb
H2(V′′)1) + D2(V′′)0) D2(V′′)1) + H2(V′′)0)
f f
VV Transfer H2(V′′)0) + D2(V′′)1) D2(V′′)0) + H2(V′′)1)
IIIa IIIb IIIc IIId IIIe IIIf IIIg
H2(V′′)1) + H H2(V′′)1) + H2(V′′)0) H2(V′′)1) + D2(V′′)0) H2(V′′)1) + He D2(V′′)1) + H2(V′′)0) D2(V′′)1) + H D2(V′′)1) + D2(V′′)0)
f f f f f f f
VT Transfer H2(V′′)0) + H H2(V′′)0) + H2(V′′)0) H2(V′′)0) + D2(V′′)0) H2(V′′)0) + He D2(V′′)0) + H2(V′′)0) D2(V′′)0) + H D2(V′′)0) + D2(V′′)0)
IVa IVb IVc
D2(V′′)1) + H D2(V′′)1) + H D2(V′′)0) + H
f f f
Reactive HD(V′′)0) + D HD(V′′)1) + D HD(V′′)0) + D
2.00(+03) 7.00(+03) 1.02(+01)
29 29 29
Va
H + H + H2
f
Recombination H2(V′′)0) + H2
2.90(+03)
30
6.00(+03)a 6.02(+01)
this work 25
3.73(+04) 7.80(+01) 1.30(+02) 1.56(+01) 1.00(+02) 1.10(+05) 8.40(+00)
24 28 26 28 27 28 26
Numbers in parentheses denote powers of 10.
one spatial direction. This leads to two tridiagonal systems of equations being solved successively with direct and efficient algorithms.34 Within each time step, variable splitting solves for the variables one by one, instead of all at one time. The updates of the first touched variables are considered by the following ones. The resulting equation system avoids the occurrence of block iteration matrices, which are difficult to handle because of their unknown properties, such as, e.g. matrix condition, diagonal dominance, eigenvalues, etc. The point iteration matrices for one variable are built up from the analytically linearized Navier-Stokes equations. These equations, which are fairly complicated to set up, are calculated by means of an algebraic programming language (MAPLE) and allow the direct calculation of the entries of the linearized Jacobi matrix.35 Compared to previous schemes, which calculated these entries by differential approximations, the present method is faster by more than a factor of 200. The coupling of the variables and the nonlinearity are taken into account within a fixed point iteration. The convergence behavior of this fixed point iteration is used to control the time step size. The degree of convergence within the fixed point iteration can be used to control the time accuracy of the integration. By tightening the tolerance parameter, the coupling iteration will be performed more thoroughly and, therefore, to a more time-accurate solution. To reach the steady-state solution as quickly as possible, the time step keeps close to the stability limit but with permanent intention to increase. Other acceleration techniques were also applied such as integration relaxation and residual smoothing.32 Integration in time can be controlled with an implicitness parameter β, which is to be set equal to 0 for explicit, 0.5 for Crank-Nicolson, and 1 for implicit integration.34 With β ) 1 the algorithm obtains the highest stability and time damping, which speeds up the approach to the steady state. The present technique permits to solve the problem as if it were fully coupled and approach the stationary solution by simultaneously avoiding complications from “stiff” flow reaction systems. The error introduced through the operator splitting technique thus can be reduced to a negligibly low value by the fixed point iteration method. The unique feature of the present
technique is that the employed defect correction method accompanying the separate iterations in each direction always reaches a stable convergence limit. This is equivalent to obtaining the fully coupled stationary solution. The code was validated through comparison with analytical solutions of simple transport diffusion reaction equations given in ref 23 and by comparison with experimental data (see Results section). The techniques presented here are specially adapted for the present problems encountered in the hydrogen system for weak reactive flow with low heat release and fast diffusion processes. In this sense it is an extreme example. The methods will not work in systems, such as combustion processes, with steep temperature gradients and fast kinetics, which cannot be decoupled from the flow field. Comparison of Experimental Data and Numerical Simulation. In the CARS measurements, due to finite spatial resolution along the beam path a radial CARS signal profile contains spatially averaged concentration information through a convolution of the square of the radial concentration profile n2(r) and a spatial instrument function f(r), which depends on the characteristics of the focused beam geometry:
Ias(p) ) ∫n2(r)f(r-p) dr
(12)
In this expression r is the radial coordinate and p the common location of the beam foci along the beam axis. The function f(r-p) quantifies the contribution to the CARS signal created at location r while the CARS focus is situated at position p. For a comparison with calculated concentration profiles, one can either deconvolve the experimental CARS data or otherwise perform a (numerical) convolution for the simulated concentration profiles to arrive at calculated CARS signal intensities. The first method is more demanding for it needs entirely recorded radial CARS profiles, which are not available for all measurement points. Therefore, in this study the second method is used. Estimation of Reaction Rate Constants by Optimization. Simulated species concentration profiles are brought into best agreement with experimental data by the variation of up to three thermal rate constants to minimize overall deviation (for example in a least-squares norm). For the numerical simulation
Numerical Simulations in Flow Reactors
J. Phys. Chem., Vol. 100, No. 22, 1996 9329
this idea leads to an optimization problem, which searches for a best fitting input parameter set. By means of an objective quality criterion, the process of parameter fitting gives unique and reproducible results for the kinetic parameters. In the case of multiple parameters to be estimated within one experiment, the use of optimization schemes offers the only way to obtain the best fit configuration, because the mutual nonlinear dependences of each variable become too complex to be adjusted by hand. The optimization algorithm either consists of a GaussNewton method with quadratic-convergence properties, or optionally a quasi-Newton scheme.36 For the Gauss-Newton method the quality criterion can be interpreted as a sum of m squares in the quality function QF:
QF(x) )
1
m
∑fi2(x) ) 1/2||f(x)||22 2 i)1
(13)
where fi(x) denotes the deviation between the measured data point i and the solution calculated numerically. The deviations fi as well as QF(x) depend on the n-dimensional kinetic parameter vector xj (j ) 1, ..., n). The Jacobian matrix J(x) of the problem contains as elements the partial derivatives of fi into xj direction. Therefore, with respect to x the gradient of QF is given by
g(x) ) ∇QF(x) ) J(x)T f(x)
(14)
and the second derivatives (curvature, or Hessian matrix) by
G(x) ) ∇g(x) ) J(x)T J(x) + Q(x)
(15)
with
Figure 5. Upper part: CARS spatial resolution function f(r-p). Measured CARS intensities by translating a 0.1 mm thick glass sheet (open circles) and a thin hydrogen flow from a nozzle (filled circles) through the common foci of pump and Stokes laser beams, respectively. The lower part of the figure shows the nozzle design for generating the spatially thin hydrogen gas flow.
m
Q(x) ) ∑fi(x)Gi(x)
(16)
i)1
The Newton equation Gp ) -g for determining the search direction p is given by
(JTk Jk + Qk)pk ) -JTk fk
(17)
with k as the iteration index of the optimization algorithm. When Qk tends to zero during the iteration procedure, eq 17 reduces to
/2||Jkpk + fk||22 f min
1
(18)
and gives the (Gauss-Newton) search direction pk for the kth iteration. Convergence properties are of second order, even though only first-order derivatives are used. By stepping pk through the parameter space (xk+l ) xk + pk), a parameter set x* will be found, which minimizes QF. Results Spatial Resolution and Convolution. For the evaluation of experimental data the CARS spatial resolution function f(rp) (instrument function), introduced in the previous section, has to be known. In practice, this parameter can be determined by translating a thin (100 µm) glass plate along the focal region of the beams and measuring the intensity of the generated nonresonant CARS signal.37 It has been argued, however, that in gas phase measurements this procedure will not reproduce the true spatial variation of the CARS signal generation. This is due to different growth mechanisms of the macroscopic CARS signal in solid and gaseous samples. For the plate the signal mainly originates from the sample boundaries whereas the contribution generated inside the thin plate is negligible.38 On
the contrary, in the gas phase the variation of refractive index is very small and the signal generated inside the medium dominates. Quantitative estimates are difficult since the effect strongly depends on the geometrical setup and beam characteristics. Taran et al.39 pointed out that the spatial extension of f(r) measured with a glass plate will be strongly underestimated. To provide the necessary input data for the modeling we investigated two approaches for the measurement of the CARS instrument function: the translation of a glass plate as well as of a thin hydrogen gas jet emerging from a 1 × 5 mm rectangular nozzle with an air coflow and probing right above the nozzle exit (as is depicted in the lower part of Figure 5). The results from both measurements are shown in the upper part of Figure 5, where the CARS intensities of both the nonresonant signal from the glass plate and the Q(1) transition in hydrogen are plotted as a function of distance from the point of maximum intensity. The longitudinal spatial resolution, defined as the distance along the common focus of the beams where 75% of the emerging signal is generated,39 turns out to be on the order of 15 and 9 mm, as measured with the nozzle and glass slide, respectively. For reasons mentioned above, the instrument function measured by means of the hydrogen jet was considered the more reliable discription of the “true” spatial resolution function. Hydrodynamics in the Laminar Flow Tube. For welldefined initial conditions in kinetic data evaluations it is desirable that a swirl free transition from inner and surrounding laminar flow profiles into a Poiseuille type flow takes place at the exit plane of the inner tube (see Figure 1). In gas dynamic terms this means streamlined flow in the mixing region so that mixing will be dominated by diffusive processes. Experimentally the ratio xoi ) Qo/Qi of outer (Qo) to inner (Qi) gas volume flow rates can be adjusted to meet these requirements. From extensive test calculations and comparisons with experi-
9330 J. Phys. Chem., Vol. 100, No. 22, 1996
Segatz et al.
Figure 6. Calculated axial velocity profile (m/s) of the total gas flow in the mixing section of the laminar flow tube reactor. The gases enter from the lower left. Volume flow ratio xoi ) 4.5 (see text), total flow 3.96 × 103 cm3/min (6.67 × 10-5 m3/s) at STP, total pressure 5.32 × 102 Pa, temperature 300 K.
mental CARS radial concentration profiles an optimum volumetric flow ratio of xoi ) 4.5 (total flow 3.96 × 103 cm3/min at STP, total pressure 533 Pa) was found by changing this parameter in the numeric modeling. Figure 6 illustrates in a three-dimensional plot the calculated axial flow velocity field in the mixing section under these conditions. For a detailed experimental investigation of convective processes especially at the end of the inner inlet tube, part of the inner helium flow was substituted by acetone, and its spatial distribution was monitored by two-dimensional laser-induced fluorescence (LIF) measurements. Because diffusional transport of acetone is significantly slower than hydrogen, the former will mainly follow the streamlines and therefore will be a tracer of convective processes. The tracer is added through the inner flow tube (the total flow is kept fixed in all experiments), while hydrogen is added through the outer side arms (2) in the laminar mixing head (see Figure 1b); Figure 7 shows its distribution in horizontal cross sections at three locations downstream from the lower end of the tube. These images demonstrate the uniform diffusive mixing and rotational symmetry of the two merging gas flows. Further downstream (corresponding to reaction times larger than 4.7 ms) signal levels, relative to dark noise and scattered light, are too low for meaningful interpretations. Likewise, Figure 8 shows the two-dimensional radial distribution of the tracer species when added through one of the inlet tubings (2) for the outer gas flow. These images demonstrate the absence of significant azimutal convective swirl that might be caused by asymmetries in the gas flows or geometric dimensions of the mixing head. The vertical development of the flow was investigated by rotating the light sheet 90° from horizontal and exciting in a plane of the flow axis. Figure 9 shows experimental intensity traces (dots) through the image plane and perpendicular to the main flow axis. It again clearly shows the gradual diffusive mixing of the tracer gas with the surrounding gas at a location starting directly below the exit of the inner tube (1) (upper part of the figure) and 42 mm downstream. The excellent agreement between measured radial concentration profiles at different axial locations and results obtained from numerical simulations (solid lines through the data points) shows that hydrodynamics and transport within the flow reactor is well described by the model. Optimization Procedure. Figure 10 shows in a projection onto the two-dimensional parameter space the path of the optimization procedure. Both the rate constant of heterogeneous relaxation of vibrationally excited hydrogen (x axis) and of vibrational energy transfer (y axis) are optimized simultaneously
Figure 7. 2D-LIF images from acetone tracer added through the central tubing in the laminar mixing head. Shown are horizontal cross sections at three distances z from the exit of the center tubing (1) in Figure 1b. In the images the LIF intensity varies from high (white pixels) to low (dark pixels).
by determination of the slope of the quality function QF in both variables. The elongated valley reflects different singularities of the Hesse matrix. Perpendicular to this valley the minimum is found within only a few iterations. After that, optimization is continued along the valley. Finally a test step is made to test for a possible saddle point of the approached minimum. Using this optimization procedure, the rate constants in the following paragraphs were derived. Vibrational Relaxation of H2(W′′)1). A visualization of the three-dimensional spatial distribution of the mass fraction of vibrationally excited hydrogen in the mixing section (with helium supplied through the central tube) is given in Figure 11. The almost instantaneous mixing after both gas flows merge, as well as the significant back-diffusion of helium from the inner tube into the supply tubes for hydrogen is clearly visible. Vibrationally excited hydrogen then is further diluted by diffusive and convective mixing and diminishes by homogeneous and heterogeneous collisional relaxation. It has to be noted from Figure 11 that to check for the effect of backdiffusion of reactive species into the center tubing the inflow border (and corresponding start of simulation) was shifted 5 mm upstream. This caused part of the inner flow tube to be “visible” in the concentration plot (sharp step at axial coordinate equal to zero). For grid points that were located in this spatial region of the inner flow tube concentration was arbitrarily set to the average value of the incoming H2(V)1) mass fraction (0.0065). It turned out, however, that this effect does not influence the CARS data, since, because diffusion is small in comparison with the high axial convection velocity, this changed the H2(V)1) concentrations by not more than 1%.
Numerical Simulations in Flow Reactors
J. Phys. Chem., Vol. 100, No. 22, 1996 9331
Figure 10. Quality function for the two-parameter optimization procedure as applied to the reactive flow of vibrationally excited hydrogen in deuterium in the flow tube with helium as carrier gas. Two-dimensional projection of QF on the optimization parameters for wall deactivation (x axis) and gas phase vibrational energy transfer rate constant (y axis). Solid line through the valley traces optimization path.
Figure 8. 2D-LIF images from acetone tracer added through one of the two sidearms (2) in the laminar mixing head in Figure 1b. Shown are horizontal cross sections at three distances z from the exit of the central tubing. In the images the LIF intensity varies from high (white pixels) to low (dark pixels).
Figure 9. Three-dimensional representation of line cuts of acetone LIF intensity distributions (points) at various distances from the inner inlet port (1) in the vertically oriented light sheet (10 mm high). Detection is through the side wall of the quartz reactor. Upper part: section from 2 to 17 mm; lower part: section from 42 to 51 mm. Solid lines: numerical simulation.
With the typical flow conditions and vibrationally excited hydrogen present (microwave discharges switched on), concentrations were measured by tuning the Stokes dye laser to
Figure 11. 3D visualization of calculated mass fraction of vibrationally excited hydrogen in the mixing section of the flow reactor (see text). Total gas pressure 5.32 × 102 Pa; p(He) 3.72 × 102 Pa; p(H2) 1.60 × 102 Pa; total gas flow 3.96 × 103 cm3/min (6.67 × 10-5 m3/s) at STP; mean flow velocity 17 m/s; microwave power 2 × 20 W; temperature 300 K.
the Q(1) transition in the first vibrationally excited state of molecular hydrogen. Hydrogen and helium are added through the outer sidearms and the central tube into the mixing head, respectively. Figure 12 shows the comparison between measured CARS intensity profiles (points with error bars) and the results of the numerical simulation (solid lines through the data points). In Figure 12a a radial CARS intensity profile at z ) 48 mm downstream from the inner flow tube exit, and in (b) and (c) the axial relative concentration profile with the center of the probe volume located on the axis of the flow tube are presented. Along this symmetry line the initial increase of the CARS signal is due to diffusion of H2(V′′)1) from the outer regions into the center of the CARS probe volume. Subsequent decay of the signal intensity can be attributed to loss of vibrationally excited species due to heterogeneous relaxation in collisions with the reactor wall (reaction Ia in Table 1) as well as homogeneous vibration-translational (VT) relaxation in collisions with hydrogen atoms, helium and molecular hydrogen (reactions IIIa,b,d in Table 1). Measured axial and radial profiles are best fitted by the simulation (solid lines in Figure 11) with a wall deactivation probability for hydrogen equal to γw ) (1.5 ( 0.3) × 10-3. The axial CARS intensity profile measured in the flow system with turbulent mixing (cf. Figure 1a) resembles data taken in the laminar system except for the mixing region (the first few centimeters). In terms of a simplified plug flow evaluation the axial coordinate z can be converted into a reaction time t )
9332 J. Phys. Chem., Vol. 100, No. 22, 1996
Segatz et al.
Figure 12. Collisional relaxation of vibrationally excited hydrogen in the flow tube: comparison of experimental H2(V′′)1) CARS spatial concentration profiles with detailed numerical simulation results. Upper left part: schematic of flow tube section from where experimental traces have been obtained; (a) radial concentration profiles 48 mm downstream from the exit of the inner gas supply tubing (1) in the laminar mixing head. The dashed and solid vertical lines mark the location of inner and outer tube walls of the reactor, respectively; (b) measured and simulated axial concentration profiles along the centerline of the flow tube; (c) same as (b) on a logarithmic scale. Conditions are as given in Figure 11.
z/Vplug, with the mean volumetric flow velocity Vplug determined by the volume flow rate. Since the time zero of the reaction is not well-defined in the turbulent mixing head, it can only be roughly estimated with the help of the plug flow velocity and an approximate distance between the end of the gas inlet tubes and the first downstream measurement location. The unsaturated CARS signal is assumed to directly correspond to the H2(V′′)1) concentration squared, assumed to be constant along the diameter of the flow tube. From half-logarithmic plots of relative H2(V′′)1) concentrations vs reaction time (neglecting data taken in the mixing region), a total first-order rate constant for vibrational deactivation of hydrogen molecules in the flow system of kplug ) 82 ( 10 s-1 can be obtained. The contribution to the pseudo-first-order rate constant from homogeneous relaxation pathways via reactions IIIa-c is calculated to be kH ) 42 s-1. From this a wall deactivation rate constant kw ) kplug - kH ) 40 ( 10 s-1 can be extracted. Within the zerodimensional approximation kw translates into a wall deactivation probability γwplug via the Kaufmann formula4
γwplug ) kwplug2r0/Vm ) (7 ( 2) × 10-4
(19)
where r0 is the radius of the flow tube and Vm the mean thermal velocity of H2(V′′)1). The obtained plug flow value γwplug is two times smaller than the one obtained by the detailed numerical simulation. Vibrational Energy Transfer of H2(W′′)1) in Collisions with D2(W′′)0). In this experiment deuterium is added through the central tube while hydrogen enters through the outer side arms (2) with the two microwave discharges switched on. Both hydrogen and deuterium were monitored in their first excited vibrational state by tuning the Stokes laser to their respective Q(1) and Q(2) transition. Figure 13b shows the axial decay of vibrationally excited hydrogen caused mainly by VV transfer in collisions with deuterium molecules and in collisions with the wall. The
Figure 13. Collisional relaxation and vibrational energy transfer of vibrationally excited hydrogen and deuterium in the flow tube: comparison of experimental H2(V′′)1) and D2(V′′)1) CARS spatial concentration profiles with detailed numerical simulation results in the flow tube. Upper left part: schematic of flow tube section from where experimental traces have been obtained; (a) radial concentration profiles of D2(V′′)1) 48 mm from the exit of the inner gas supply tubing (1) in the laminar mixing head. Shown as dashed and solid vertical lines are the geometric dimensions of inner and outer tube walls of the reactor, respectively; (b) measured and simulated axial concentration profiles of H2(V′′)1) along the centerline of the flow tube; (c) same as (b) but for D2(V′′)1).
buildup of vibrationally excited deuterium due to VV processes and its subsequent decay mainly due to wall loss can be followed along the flow axis in part c and its radial distribution (48 mm downstream from the mixing section) in part a of the figure. The relevant reaction mechanism used in the simulations is listed in Table 1. The previously obtained value for the wall deactivation probability γw ) (1.49 ( 0.3) × 10-3 is taken for both isotopic hydrogen species,15 since within the accuracy of the measurements reported there, an isotope effect was absent. This can partly be explained by the far removed optical phonon frequencies of the fully fluorinated carbon bonds used in the type of Teflon wall coating40 from the hydrogen vibrational frequencies. Data points are best fitted from the detailed numerical simulation with the rate constant for vibrational energy transfer from hydrogen to deuterium equal to kvv ) (6 ( 0.5) × 109 cm3 mol-1 s-1. For comparison with the plug flow evaluation, again the axial coordinate is replaced by reaction time t ) z/Vplug. In this reactor data from the turbulent mixing region are not available, so the initial buildup of D2(V′′)1) cannot be observed. An attempt was made to fit concentration-time histories by results from a numerical integration of the homogeneous reaction mechanism shown in Table 1, by varying the rate constant for the vibrational energy transfer kvv. Heterogeneous relaxation was taken into account in terms of the pseudo homogeneous rate constant kwplug determined previously. The nearly exponential decay of H2(V′′)1) is best fitted by a rate constant kvv ) 4.3 × 109 cm3 mol-1 s-1, almost 40% smaller than the value obtained by the detailed numerical simulation. The D2(V′′)1) concentration profile can only be fitted when data points are shifted by 2 ms on the time axis. This means that the mean residence time of reactants in the mixing region is longer than the time estimated by the plug flow velocity and the distance between the end of the inlet tubes (Figure 1a) and the first measurement location. This mismatch can be reasonably explained by the increased
Numerical Simulations in Flow Reactors
J. Phys. Chem., Vol. 100, No. 22, 1996 9333
residence times in the mixing region caused by the turbulent motion, and again shows the importance of a detailed knowledge of the mutual dependence between reaction and flow field characteristics.
ches Rechnen in Mathematik und Naturwissenschaften” of the Interdisziplina¨res Zentrum fu¨r Wissenschaftliches Rechnen, Heidelberg.
Discussion and Conclusions
(1) Wood, R. W. Proc. R. Soc. (London) 1922, A102, 1. (2) Bonhoeffer, K. F. Z. Phys. Chem. 1924, 133, 199. Boehm, R.; Bonhoeffer, K. F. Z. Phys. Chem. 1926, 119, 385. Bonhoeffer, K. F.; Harteck, P. Z. Phys. Chem. 1928, A139, 64. (3) Schlichting, H. Boundary layer theory; G. Braun: Karlsruhe, 1982. (4) Kaufmann, F. Progress Reaction Kinetics; Pergamon Press: New York, 1961; Vol. 1. (5) Poirier, R. V.; Carr, R. C., Jr. J. Phys. Chem. 1971, 75, 1593. (6) Ogren, P. O. J. Phys. Chem. 1975, 79, 1749. (7) Westbrook, C. K.; Creighton, J.; Lund, C.; Dreyer, F. L. J. Phys. Chem. 1977, 81, 2542. (8) Keyser, L. F. J. Phys. Chem. 1984, 88, 4750. (9) Walker, R. E. Phys. Fluids 1961, 4, 1211. (10) Brown, R. L. Res. Natl. Bur. Stand. (U.S.) 1978, 83, 1. (11) Abbatt, J. P. D.; Demerjian, K. L.; Anderson, J. G. J. Phys. Chem. 1990, 94, 4566. (12) Seeley, J. V.; Jayne, J. T.; Molina, M. J. Int. J. Chem. Kinet. 1993, 25, 571. (13) Billing, G. D.; Kolesnick, R. E. Chem. Phys. Lett. 1993, 215, 571. (14) Buchenau, H.; Toennies, J. P.; Arnold, J.; Wolfrum, J. Ber. Bunsenges. Phys. Chem. 1990, 94, 1231. (15) Arnold, J.; Bouche´, T.; Dreier, T.; Wichmann, J.; Wolfrum, J. Chem. Phys. Lett. 1993, 203, 283. (16) Seitzman, J. M.; Hanson, R. K. In Instrumentation for Flows with Combustion; Taylor, A. M. K. P., Ed.; Academic Press: London, 1993. (17) Lozano, A.; Yip, B.; Hanson, R. K. Exp. Fluids 1992, 13, 369. (18) Rahn, L. A.; Zych, L. J.; Mattern, P. L. Opt. Commun. 1979, 39, 249. (19) Clyne, M. A. A.; Thrush, B. A. Proc. R. Soc. 1963, A275, 544. (20) Hartley, D. B.; Thrush, B. A. Proc. R. Soc. 1967, A297, 520. (21) Pe´alat, M.; Levebvre, M.; Taran, J. P. E.; Kelley, P. L. Phys. ReV. A 1988, 38, 1948. (22) Oran, E. S.; Boris, J. P. Numerical Simulation of ReactiVe Flows; Elsevier Science Publ. Inc.: New York, 1987. (23) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; J. Wiley & Sons Inc.: New York, 1960. (24) Schatz, G. C. Chem. Phys. Lett. 1984, 108, 532. (25) Arnold, J.; Dreier, T.; Chandler, D. W. Chem. Phys. 1989, 133, 123. (26) Lukasik, J.; Ducuing, J. Chem. Phys. Lett. 1974, 27, 203. (27) Calvert, J. B. J. Chem. Phys. 1972, 56, 5071. (28) Bowes, C.; Mina, N.; Teitelbaum, H. J. Chem. Soc., Faraday Trans. 1991, 87, 229. (29) Sun, Q.; Bowman, J. M. J. Chem. Phys. 1990, 94, 718. (30) Trainor, D. W.; Ham, D. O.; Kaufmann, F. J. Chem. Phys. 1973, 58, 4599. (31) Gershenzon, Yu. M.; Ivanov, A. V.; Rozenshtein, V. B.; Umanskii, S. Ya. Prog. React. Kinet. 1961, 1, 1. (32) Anderson, D. A.; Fannehill, J. C.; Pleteher, R. H. Computational Fluid Mechanics and Heat Transfer; Series in Computational Methods in Mechanics and Thermal Sciences; Bristol, PA, 1984. (33) Issa, R. I. J. Comput. Phys. 1985, 62, 40. (34) Fletcher, C. A. J. Computational Techniques for Fluid Dynamics; Springer: Heidelberg, FRG, 1988; Vol. 1, 2. (35) Segatz, J. Dissertation, University of Heidelberg, 1995. (36) Gill, P. E.; Murray, W.; Wright, M. H. Practical Optimization; Academic Press: London, 1992. (37) Boquillon, J. P.; Pealat, M.; Bouchardy, P.; Collin, G.; Magre, P. Opt. Lett. 1988, 13, 722. (38) Bloembergen, N. Nonlinear Optics; W. A. Benjamin, Inc.: New York, 1965. (39) Druet, S. A. J.; Taran, J. P. E. Prog. Quant. Electron. 1981, 7, 1. (40) Gershenzon, Yu. M.; Ivanov, A. V.; Kucheryavi, S. I.; Lyapunov, A. Ya.; Rozenshtein, V. B. Kinet. Catal. 1986, 27, 928. (41) Pirkle, R. J.; Cool, T. A. Chem. Phys. Lett. 1976, 42, 58. (42) Bott, J. F. J. Chem. Phys. 1976, 65, 3921. (43) Cacciatore, M.; Billing, G. D. J. Phys. Chem. 1992, 96, 217.
Comparing the plug flow results with those obtained by a detailed numerical simulation within a suitable reactor, it is evident that, even at low pressures within the plug flow limit, the sum of neglected aspects can lead to significant underestimations of evaluated rate constants, approaching the 40% range. The Kaufmann equation (19) is derived under the assumption that the number of molecules, ∆N, deactivated at the surface S ) 2πr∆z in a flow tube with radius r and length element ∆z within time ∆t, ∆N ) (1/4)VmnSγw∆t, is distributed homogeneously within the volume V ) πr2∆z, giving rise to a uniform decrease in concentration ∆n ) ∆N/V ) Vm/(2r0) γwn∆t without radial concentration gradients. This holds for small γw and fast diffusive transport. However, it is evident from Figure 9 or 11 that the concentration profiles do show radial variations giving rise to a systematic error in derived rate constants, if the deactivation probability is evaluated in terms of a homogeneous process. Due to the nonuniform radial concentration distribution (concentrations being higher in the high-speed domains in the center of the reactor and lower in the low-speed domains near the wall) the average axial transport of excited species takes place with an effective average transport velocity larger than the plug flow velocity. This effect is further enhanced in our measurements of axial concentration profiles because CARS detection is more weighted towards the location of the beam foci (i.e., on the centerline of the reactor), so that molecules in high-speed domains mainly contribute to the total CARS signal. In an assumed “worst case” of a parabolic concentration profile as a result of less effective diffusion, and detection limited to one point located on the reactor axis, the measured effective average transport velocity will equal the maximum of the Poiseuille flow velocity which, in this case, is twice the plug flow velocity. Since velocity directly enters into the evaluated time scale, rate constants determined by plug flow evaluations are systematically lower. The rate constant for the VV transfer presented in this work is in good agreement with the results of Pirkle and Cool41 and Bott42 of 5.9 × 109 and 6.1 × 109 cm3 mol-1 s-1, respectively. They employed vibrational excitation of hydrogen via VV transfer in collisions with vibrationally excited HF. The subsequent relaxation of H2(V′′)1) due to VV transfer to D2(V′′)0) is monitored via time resolved HF fluorescence. The experimental results for kvv can be opposed with the most recent semiclassical calculations of Cacciatore et al.43 and Billing et al.:13 They obtained values of 2 × 109, 2.8 × 109, and 3.1 × 109 cm3 mol-1 s-1 from SCF and CI calculations, respectively. This shows that a quantitative ab initio modelling of even the simplest VV energy transfer is still lacking. Acknowledgment. This work is supported by the Deutsche Forschungsgemeinschaft (SFB 359). J.S. is thankful to his grant from the Graduiertenkolleg “Modellierung und wissenschaftli-
References and Notes
JP953578H
