9794
Ind. Eng. Chem. Res. 2005, 44, 9794-9804
Modeling of Multicomponent Concentration Profiles in Membrane Microreactors Khaled Alfadhel and Mayuresh V. Kothare* Integrated Microchemical Systems Laboratory, Department of Chemical Engineering, Iacocca Hall, 111 Research Drive, Lehigh University, Bethlehem, Pennsylvania 18015
In this paper, we study the problem of modeling multicomponent concentration profiles in a membrane-based microreactor. Using basic constitutive laws of mass balance, we derive a loworder model of a generic membrane microreactor, incorporating chemical reaction and permeation through a selective micro-membrane and utilizing a pressure distribution formula for slipping flows from our previous work. Without loss of generality, the model can address nonisothermal conditions and can be extended to allow flow compressibility. We study the utility of our model in evaluating the optimal design and operation of a palladium-based membrane microreactor for conducting hydrogen purification and water gas shift (WGS) reforming in microfuel processing applications. 1. Introduction Microchemical systems are a new generation of miniature chemical systems that perform chemical reactions and separations in precisely fabricated threedimensional microreactor configurations in the size range of a few micrometers to a few hundred micrometers.1 Among the distinguishing features of microchemical systems are their high surface area per unit volume and small diffusion paths, both of which provide opportunities for achieving performance characteristics that are significantly different from classical full-scale chemical processes. Once standardized, these automated individual microchemical systems can be replicated in parallel, leading to a microchemical plant “flow sheet” with individually optimized identical micro-units.1 The focus of this paper is on modeling membranebased microreactors, which are an important class of microreactors that can combine reaction and separation in one single microdevice.2-5 Applications of membrane microreactors include the dehydrogenation of cyclohexane to benzene,3 hydrogen generation by oxidative coupling of methane,4 the production of moisture-free formaldehyde by the dehydrogenation of methanol,6 and hydrogen separation/water gas shift (WGS) reaction.2,5 Prototype membrane-based microreactors,2,5 as well as numerous other microchemical systems, have been fabricated and are now being used in a wide range of chemical and biological applications.7 However, despite numerous advances, the fabrication of microchemical systems remains extraordinarily difficult, nonstandard, and based on customized designs that use microfabrication techniques borrowed from the microelectronics industry. All of these factors make reproducibility of results and experimental validation of microreactor performance quite tedious, if not impossible. Therefore, there is a growing need for proper mathematical modeling that can describe the physical phenomena and operability of these systems with ease. Adequate models, which involve many physical parameters that can only * To whom correspondence should be addressed. Tel.: (610) 758-6654. Fax: (610) 758-5057. E-mail:
[email protected].
be found experimentally, are classically complicated because of the coupling of heat, mass, and momentum transfer in confined geometries. Therefore, the theory and fundamental physical understanding of microchemical systems remains an area of ongoing research.1,8 The goal of this paper is to develop a low-order mathematical model to describe multicomponent concentration profiles in a membrane-based catalytic microreactor/microseparator and to evaluate its utility in predicting performance of a WGS microreactor as part of an integrated micro-fuel-cell power device. Our objective is to provide a generic model that includes microscale phenomena while maintaining reasonable mathematical solvability for possible use in optimization and design studies. Simulations based on computational fluid dynamics have been utilized in the literature to characterize results found from specific experiments involving specific microgeometries. Thus, it is difficult to use these models for different geometries, other transport phenomena, or to optimize the performance of other microchemical systems.9 On the other hand, molecular-based approaches are known to be computationally very expensive and require a detailed understanding of surface, fluid, and surface-fluid characteristics,10 which is often not easily available. Alternatively, one-dimensional models that address the important phenomena in microscale processess have been shown to be sufficient in providing reliable results.11 We adopt this philosophy in our modeling approach, based on the conclusions presented in refs 9 and 11. The mathematical model we derive describes multicomponent concentration profiles in microchannels with rectangular cross-sections in which both a nonelementary, non-constant-volume, reversible heterogeneous chemical reaction and selective membrane permeation occur. The model incorporates the pressure drop formula that we derived in our previous work12 for slipping flows in rectangular microchannels, and the model is made general to treat the permeating species as a reactant, a product, or even an inert. Without loss of generality, the model can handle isothermal and nonisothermal conditions.
10.1021/ie049176q CCC: $30.25 © 2005 American Chemical Society Published on Web 11/03/2005
Ind. Eng. Chem. Res., Vol. 44, No. 26, 2005 9795
2. Literature Review Many attempts on the modeling of macroscale catalytic membrane reactors have been reported in the literature.13-19 Hara et al.13 developed an ideal gas model that incorporates mixing diffusion in the radial direction of a catalytic packed-bed reactor. Kikuchi et al.14 formulated a model that excluded radial concentration gradients in a packed-bed membrane reactor. The effect of the sweep gas on permeation was modeled by Barbieri and Di Maio.15 Three different steady-state models for traditional fixed-bed reactor, a catalytic mesoporous ceramic membrane reactor, and a palladium membrane reactor were developed and compared by Criscuoli et al.16 Barbieri et al.17 studied methane reforming in a palladium-based membrane reactor and concluded that conversion increased as the temperature, feed flow rate, and steam-to-methane molar ratio each increased. The dehydrogenation of cyclohexane to benzene in catalytic membrane reactors was performed by Ito et al.18 and Moon et al.19 Component concentration profiles were computed by Ito et al.,18 based on a differential shell model under isobaric conditions. Along with the flow model, Moon et al.19 incorporated a pseudo-homogeneous model for the catalytic bed. All of the previous models assumed plug flow, isothermal and isobaric conditions, and were derived for macroscale applications. Previous work reported in the area of microscale catalytic membrane reactors is largely experimental.2,5,20 Franz et al.2 integrated a palladium-based membrane microreactor with heaters and sensors. Karnik et al.5,20 demonstrated H2 permeation in an ultrathin palladium micromembrane. Experimental membrane-based microreactor prototypes were also reported in refs 3 and 4. To the best of our knowledge, no effort has been made to develop mathematical models for microscale membrane reactors. Microscale membrane reactors can, in principle, be modeled using continuum and noncontinuum approaches (such as equilibrium and nonequilibrium molecular dynamics21,22) or kinetic techniques (such as the direct simulation Monte Carlo method23-27). However, molecular dynamic approaches, despite their ability to provide results that are in excellent agreement with experimental measurements, are applicable only for very small liquid volumes with dimensions and time intervals on the order of nanometers and nanoseconds, respectively.10 Kinetic techniques, on the other hand, have a wider range of validity21,22 but (i) are inevitably computationally expensive26 and (ii) are known to have numerous limitations and disadvantages.10 An issue such as the thermal velocity in microchemical systems is capable of creating significant noise in the simulation. This noise requires extremely large sample sizes to obtain physically reliable results, even for the simplest flow problems.26 To overcome the many difficulties in formulating molecular models and solving them, temporal and/or spatial quantities are averaged to construct solvable models.23 Simplified numerical approaches and analytically derived expressions were reported to be capable of reproducing molecular simulation predictions in microscale processes.24-27 In the work of Lı´sal et al.,24 the reaction-ensemble Monte Carlo simulation was observed to agree with two different conventional thermodynamic approaches. Coppens and Malek25 reported an excellent similarity when they compared their Monte Carlo
Figure 1. Schematic representation of the microchannel to be modeled.
results to the relevant analytical expressions derived under the continuum assumption for reactions in nanopores. Snyder et al.27 demonstrated that gradient continuous time Monte Carlo simulations were in excellent agreement with their mesoscopic solutions for microporous membranes. Shen et al.26 compared molecular approaches to the solutions of Navier-Stokes equations for rarefied gas flows in microchannels and verified that, for Knudsen numbers of Kn < 0.1, continuum solutions were valid. To model microscale phenomena using continuum approaches, appropriate adjustments to the boundary conditions must be incorporated.10 Because of the proven applicability of continuum models for Kn < 0.1, which is the limit most microchemical systems do not exceed, our previous work12 evaluated the impact of wall velocity slip and selective component permeation through the membrane on the velocity profile. Our model solutions were derived in an analytical form, using finite Fourier transforms, and were shown to recover several known results from the literature as special cases. In light of their powerful utility in providing reliable results,10 we use continuum principles to model the composition profiles in microscale membrane reactors. 3. Mathematical Model Multicomponent mass-transfer problems are known to be inherently complex and difficult to model. There is no single model in the existing literature that fully describes and addresses the wide spectrum of multicomponent modeling problems. Only a few exact or approximate solutions to some special situations exist in the literature. For the gas phase, the theory of Maxwell-Stephan is believed to be a rigorous starting point.28 Maxwell-Stephan equations are used when diffusion is the only driving force.29 However, there are other applications in which fluid flow as well as temperature gradients become a means for mass transfer. These applications require more complicated models, which are described by partial differential equations in three dimensions (3D). In microscale systems, the larger the number of boundary conditions in a transport problem, the more complex the problem becomes. This is due to the modification that each boundary condition must undergo to address the secondary effects associated with microscales. In this paper, we adopt a mathematical approach to model the multicomponent concentration profiles in membrane microreactors based on a simpli-
9796
Ind. Eng. Chem. Res., Vol. 44, No. 26, 2005
fication that circumvents this problem and requires no modification of the boundary conditions. With this approach, we provide a computationally tractable alternative to molecular modeling formulations that gives predictive tools that are easier to apply. This approach is also well-justified based on the recent results reported in refs 9 and 11. Figure 1 shows a schematic of the prototypical membrane microreactor under consideration. A gaseous mixture enters the microchannel on the left and flows axially. The top wall of the microchannel is assumed to be a membrane that is permeable to one species in the mixture. Either the top or bottom wall of the microchannel is assumed to be catalytically activated by deposition or coating of a catalyst layer for the heterogeneous reaction to occur. We start with the individual species balance equation in the membrane microreactor:
rate of accumulation ) in - out ( reaction permeation (1) Equation 1 is valid at each point; however, what is important here is the formulation of a continuous model that gives the concentration trajectory relating the mole fraction at any point to a known reference such as the inlet values. To reach that goal, we state the definition of xi under steady-state conditions:
xi )
(Jin i )Acs + ξi
∫A
r
∫A
Jr dAr - qi
r
m
p
Jp dAp + ∆ξ
r
(3)
Because of the rectangular cross-sectional area of the microchannel, uniform membrane coverage, and the even catalyst exposure, there will be no concentration gradient in the x-direction. Because permeation through the membrane is at least 1 order of magnitude less than diffusion through the fluid bulk, permeation can be viewed as the rate-determining step. Thus, gradients in the z-direction can be ignored.14,15 Incorporating mixing diffusion in the radial direction was not found to greatly reduce the discrepancy between model outcomes and the experimental findings of Hara et al.13
∫
∫
Equation 4 gives the cross-sectional average mole fraction at an axial position y in an integral form. For most applications, the differential form is preferred. To derive the differential version, we perform the cross multiplication and differentiate twice, with respect to the spatial dimension y. After doing some algebraic manipulations, the final result reduces to
( ) ( ) [
0 ) 2[(∆ξ)Jr - Jp] d2xi dy2
dxi dy
2
+ [(ξi - (∆ξ)xi)Jr + (xi -
]( )
dJr dJp dxi + (xi - qi) dy dy dy (5)
- (ξi - (∆ξ)xi)
(2)
∫A Jr dAr
∫
y y T [q J dy - ξi 0 Jr dy] Pνave in i 0 p y y R T 1[ J dy - ∆ξ 0 Jr dy] Lz Pνave in 0 p (4)
xi(y) )
ApJp dAp
Qm
Jin ∑ i ) - ∫A i)1
( ) ( )∫
xi,in - (R/Lz)
qi)Jp]
where Jin i is the inlet flux of species i and qi is an indicator that assumes a value of unity if species i permeates and zero otherwise. Jp and Jr are the permeation and reaction molar fluxes, respectively, and Qm is the total molar flow rate. As shown in Figure 1, Acs ) 2LxLz, Ar ) (2Lx)y, and Ap ) (2Lx)y are the crosssectional area of the rectangular microchannel, the reaction area, and permeation area, respectively, in terms of the microchannel width (2Lx), depth (Lz), and the spatial length (y). Depending on the stoichiometric coefficients ξi of the reactants and products, the reaction volume may change, depending on the sign and the value of ∆ξ. (Note that ∆ξ is the difference between the sum of the stoichiometric coefficients of the products and that of the reactants. Thus, ∆ξ ) ∑prodξi - ∑reactξi.]. The permeation flux will constantly decrease the overall molar flow rate. Thus, the total molar flow rate can be written as
Qm ) Acs(
Vertical gradients can be critical when slow diffusion problems, two-phase flows, or fast reactions are involved. With the ideal gas assumption,13-17 and utilizing the fact that mole fractions sum to unity, eq 2 can be written in terms of the species inlet mole fraction (xi,in), the average stream velocity (υave), the absolute temperature and pressure (T and P, respectively), and the universal gas constant (R):
Equation 5 is a second-order ordinary differential equation (ODE) that was derived from an integral model. Because of differentiation, the solution generated by this differential model loses its uniqueness. Uniqueness can be ensured by defining a proper set of boundary conditions to recover the integral model. The only set of boundary conditions available is the set of inlet mole fractions. Without an additional set of boundary conditions, the problem remains underspecified. Imposing improper boundary conditions, such as no gradient at the inlet or constant concentration at an infinite length of the microchannel, will cause model singularity and will consequently lead to solution discontinuities, as will be shown later in this work. The first improper set of boundary conditions makes no sense physically, whereas the second set is computationally very expensive. Furthermore, when one applies a boundary condition at the channel outlet, the model becomes a boundary-value problem, which is more difficult to solve. Therefore, the proper set of boundary conditions must be addressed differently. Starting from eq 1, one can write the following equation for any species:
(
) (
)
Pνavexi Pνavexi ) T T
in
)
R [ξ Lz i
∫0y Jr dy - qi∫0y Jp dy] (6)
and, from eq 3, with the ideal gas assumption, one can write the following ODE:
( )
R d Pνave ) [(∆ξ)Jr - Jp] dy T Lz Differentiating eq 6 and substituting eq 7, we get
(7)
Ind. Eng. Chem. Res., Vol. 44, No. 26, 2005 9797
(
)
dxi RT [(ξi - (∆ξ)xi)Jr + (xi - qi)Jp] ) dy LzPνave
(8)
Equation 8 is now the second boundary condition needed for the well-posedness of the problem. Specifying T, P, υave, and the inlet mole fractions, and given the reaction and permeation rate expressions, the derivative is defined at the inlet of the microchannel. The model is now an initial value problem that can be solved by numerous integration techniques. To show the advantages of using microreactors for gaseous flow applications, eq 8 can be written in terms of the hydraulic diameter (hD), the molecular weight of the stream (MW), and the average viscosity (µ):
(
)
dxi 1 MW/µ ) [(ξi - (∆ξ)xi)Jr + (xi - qi)Jp] (9) dy Re Lz/hD With the Reynolds number (Re) being in the denominator, the microreactors clearly exhibit larger gradients for small Re and, thus, guarantee compactness. Typical Re values for microdevices are well below 100, as opposed to large-scale applications, which normally operate at much higher Re and even in turbulent flows. For low-Re flows, it is important to notice that the no-gradient boundary condition is inappropriate for microchemical applications, because the first derivative at the inlet vanishes only when
() Jr Jp
)
in
qi - xi,in
(10)
ξi - (∆ξ)xi,in
In general form, the first derivative vanishes and the solution experiences a mathematical singularity whenever xi satisfies the following equation:
xi )
(11)
Jp - (∆ξ)Jr
3.1. Model Reduction. Before we show the numerical approaches with more general forms of the reaction and permeation fluxes, it is noteworthy to list the reduced-order model equations for the two extreme cases of no reaction (Jr ) 0) and no permeation (Jp ) 0). When Jr ) 0 and Jp * 0, eq 5 becomes
[ ( )]
)
d [ln(Jp)] dy (12)
Upon integration once, using eq 9 as a boundary condition, one can get the first derivative to be
( ) [
)]
dxi hD MW/µ (xi - qi)2 ) J dy Lz Re in xi,in - qi p
(
[
xi,in - qi
)]
hD MW/µ Lz Re in
(
∫xx
i
dθ (θ - qi)2Jp
i,in
xi,in
∫0yJp dy
1 - xi,in{(hD/Lz)[(MW/µ)/Re]in}
(14)
For the nonpermeating species (qi ) 0), the species mole fractions can be evaluated as follows (note that
(15)
Generally, equating eq 9, which is valid at any point in the micromembrane separator, to eq 13 leads to the following result:
xi ) qi +
MW/(µRe) (x - qi) MW/(µRe)in i,in
(16)
Similarly, when Jp ) 0 and Jr * 0, eq 5 becomes
( ) [
)]
dxi hD MW/µ [ξi - (∆ξ)xi]2 ) J dy Lz Re in ξi - (∆ξ)xi,in r
(
(17)
which leads to a similar integral solution:
y)
ξi - (∆ξ)xi,in
∫xx
dθ (18) [ξi - (∆ξ)θ]2Jr
i
(hD/Lz)[MW/(µRe)in]
i,in
Given the observation of the similarity between eqs 9 and 17, one can find the following result for ∆ξ * 0:
xi )
(
MW/(µRe) 1 ξ [ξ - (∆ξ)xi,in] ∆ξ i [MW/(µRe)]in i
)
(19)
When ∆ξ ) 0, the species mole fraction in a microreactor can be evaluated from the following equation:
[(
) ]∫
hD (MW/µ) Lz Re
in
y
0
Jr dy
(20)
A model is not needed for the case when Jp ) Jr ) 0, because the concentrations of all species will remain invariant, under the assumption of negligible thermal diffusion. 3.2. Illustrative Case. The solution of the secondorder model gives the cross-sectional average concentration profile of any species, as a function of temperature, total pressure, and the other species’ concentrations. To illustrate this result, we will assume an ideal gas mixture of five components: A, B, C, D, and I. The first four components undergo the following reversible nonelementary heterogeneous reaction that occurs at z ) 0 or z ) Lz, the lower or upper surfaces of the microchannel, respectively: k1
(13)
which has an implicit general analytical solution in the form of the following integral equation:
y)
xi )
xi ) xi,in + ξi
qiJp - iJr
dxi d d 0 ) 2 [ln(xi - qi)] ln dy dy dy
the maximum value of the term ({(1/Re)[(MW/µ)/(Lz/ Dh)]}in ∫0y Jp dy) is the inlet mole fraction of the permeating species):
z ξCC + ξDD ξAA + ξBB y\ k
(21)
2
with ∆ξ ) (ξC + ξD) - (ξA + ξB). Assuming that the Arrhenius-type reaction rate constants k1 and k2 account for all heterogeneous effects,15 the ideal gas reaction rate expression, in terms of the species concentration Ci, will be
Jr ) r ) k1CRACβB - k2CγCCδD ) k1xRAxβB
R+β
(RTP )
- k2xγCxδD
(RTP )
γ+δ
(22)
9798
Ind. Eng. Chem. Res., Vol. 44, No. 26, 2005
To allow for nonelementary reactions, the exponents R, β, γ, and δ of the concentrations in the rate expressions will be allowed to be different from the stoichiometric coefficients. Component I is an inert, from a reaction perspective. We further assume that the membrane is permselective to one component, S, which can be only one of the five species. Because of the permselectivity, the flag qi was introduced to have a value of unity for the permeating species and zero otherwise. Membranes in microsystems are known to be ultrathin, for which permeation is typically expressed as follows:
Jp )
k hp n [P - P h n] φ s,2 k hp )
where φ is the membrane thickness, k h p the permeability coefficient, and kp the Arrhenius-type permeation coefficient (with a permeation pre-exponential factor kp,0 and an activation energy Ep). Ps,2 is the partial pressure of the permeating species on the high-pressure side (the gaseous mixture), and P h is its partial pressure on the low-pressure side (the sweep gas). The parameter ξ was introduced to allow more flexibility in the permeation mechanism. If diffusion through the membrane is the rate-controlling step in the permeation mechanism, the permeation coefficient will be purely Arrhenius type ( ) 0). However, there are some permeation mechanisms that are controlled by other steps for which permeation coefficients might follow other non-Arrhenius-type expressions. One example is the palladium membrane used in our previous work, which required a value of ) 1.12,30 Depending on the thickness of the membrane and the permeation mechanism, the exponent n can have values from n ) 0.5 (Sievert’s law) to n ) 1 (Fickian diffusion). For these general permeation and reversible reaction expressions, both dJp/dy and dJr/dy must be evaluated to arrive at the final derivative form of eq 5. These derivatives require numerous algebraic steps before they can be substituted in the model, because of the generality of allowing nonisothermal and nonisobaric slip conditions, as well as variable partial pressure on the low-pressure side of the membrane. Letting E1 and E2 represent the forward and reversible activation energies, respectively, xs denote the mole fraction of the permeating species, and K h ) (P h /xsP)n, permeation occurs when xsP > P h , thus, kp/φT (xsP)n > 0. Dividing by the latter term to define new sets of concentration-dependent dimensionless numbers (Damko¨hler numbers Da1, Da2, and Da) and differentiation of eqs 22 and 23 leads to the final second-order model: 2
+
h )] [(ξi - (∆ξ)xi)Da + (xi - qi)(1 - K
( ) ( ) d2xi dy2
F(y, T, T′, P, P′, P h, P h ′, xi, x′i, xj*i,x′j*i) where
[kp/(φT)](xsP)n
) Da1
xRA xβB xsn
k2xγCxδD[P/(RT)]γ+δ xγC xδD ) Da2 n kp xs (xsP)n φT Da ) Da1 - Da2
(26)
(27)
(28)
and the function F(y, T, T′, P, P′, P h, P h ′, xi, x′i, xj*i, x′j*i) is given as
(24)
T
( )
Da2 )
k1xRAxβB[P/(RT)]R+β
(23)
kp
dxi 0 ) 2[(∆ξ)Da - (1 - K h )] dy
Da1 )
-
dxi (25) dy
with the boundary conditions
( ) dxi dy
)
y)0
xi(0) ) xi,in
[
(30)
)]
hD kp(xsP)n MW/µ [(ξ - (∆ξ)xi,in)Dain + φLz Re in i T h in)] (31) (xi,in - qi)(1 - K
(
3.2.1. Asymptotic Models. The cases when either or both Da1 and Da2 approach 0 or ∞ are important from both physical and mathematical standpoints. The first extreme is when Da1 ) Da2 ) 0, which corresponds to the case when the rate of the heterogeneous reaction is insignificant, with respect to permeation (Jr f 0). Setting Da1 ) Da2 ) 0 does not alter the order of the differential equation; thus, regular perturbation solutions are valid. For Da1 ) Da2 ) 0, the permeation model becomes
0)
[
( )
dxi (1 - K h) 2 n xi - qi dy
(
-
)
2
-
( )] { d2xi dy
2
+
d ln(xsP) dy
}( )
Ep (1 - K h ) d ln T d ln P h dxi (32) -K h RT n dy dy dy
Because dxi/dy * 0, the previously derived eq 13 can be solved as a reduced-order model with the Jp expression from eq 23. (Note: When dxi/dy ) 0, xi can be evaluated as described during the model derivation.) Equation 32 represents a system of nonlinear ODEs and can, with some algebraic manipulations, be reduced to the previously derived model given as eq 13. Any analytical approach should start with species S before addressing the equations for the other species. No discontinuity is expected, for the same reason previously discussed. Because of the insignificance of the heterogeneous reaction, the terms that came from the reaction rate expression disappeared. Having solved for S, the other species concentration profiles can be found from eq 15. Similarly, when permeation is insignificant, compared to either the forward (Da1 f ∞) or the reverse (Da2 f ∞) reactions, eq 17 will be the reduced order form of the
Ind. Eng. Chem. Res., Vol. 44, No. 26, 2005 9799
model with Jr being the faster rate for which Da approaches ∞. The advantage of solving reduced-order models is that they only require inlet concentrations to solve the systems of first-order differential equations. (Note that, without any loss of generality, the model can be extended to any number of species.) These equations have similar form that makes it possible to derive the concentration profiles of the other species once the solution to one component is known. For Jr ) 0 and from eqs 13, one can obtain the following:
( )
dxi xj,in - qj xi - qi ) dxj xi,in - qi xj - qj
2
(33)
Solving eq 33, with xi(xj,in) ) xi,in, gives a solution of this form:
xj ) qj +
xj,in - qj (x - qi) xi,in - qi i
(34)
Because one must solve for the permeating species first, the rest of the species are evaluated without the need for solving a system of differential equations:
xj ) xj,in
(
1 - xs 1 - xs,in
)
(35)
Equation 17 can be handled in a similar fashion. When Jp ) 0 and for ∆ξ ) 0,
ξj xj ) xj,in - (xi,in - xi) ξi and when ∆ξ * 0,
xj )
( [
]
ξj - (∆ξ)xj,in 1 ξj (ξ - (∆ξ)xi) ∆ξ ξi - (∆ξ)xi,in i
(36)
)
(37)
When solving for an inert, eq 37 gives
xI ) xI,in
(
ξi - (∆ξ)xi
ξi - (∆ξ)xi,in
)
(38)
which predicts that xI ) xI,in when ∆ξ ) 0. This is not unexpected, because this model dictates that permeation is insignificant, compared to the chemical reaction, and, because the inert species is not part of that reaction expression, i.e., ξI ) 0, its concentration shall not be affected under the assumption of constant reaction volume. The same conclusion can be reached when eq 36 is solved for an inert. It is noteworthy that the solutions to eqs 34, 36, and 37 are not only independent of the rate expressions, but are also linear. One mole fraction besides the inlet mole fraction of species i fixes the overall concentration profiles of species j. Because of its inherent degree of nonlinearity, the general reactive/separation model (eq 25) should be discretized for the numerical solution to be generated. The pressure drop formula we determined in our previous work12 can be substituted in the model along with our permeation results.5,12 Since we have not yet considered any temperature effects, the model will be solved under isothermal conditions, although it was derived for the general case of nonisothermal conditions. The model represents a set of four independent second-
order nonlinear ODEs. The last species can be found from the sum of the mole fractions. The most efficient way to solve the model is to generate the equations for reacting species (A, B, C, and D), coupled with the pressure and temperature drop equations and then evaluate the solution. Upon convergence, the profile of the last species is found using the sum of the mole fractions. 4. Case Study: Palladium Membrane Microreactor Because of the nonlinearity of our model, numerical solutions are typically difficult.10 Taking advantage of the theoretical observations we made throughout the derivation of the model, lower-order models can be found. From the inlet boundary conditions, one can notice that eq 31 was not limited to the inlet of the channel, and thus can be used as a first-order model for the concentration profiles instead of eq 25. However, there is an added expense of having the velocity function, which is an integral equation that must be evaluated at each mesh point after the mole fraction of the permeating species has been determined. Thus, the first-order model can be written as
( ) [
)]
hD kp(xsP)n MW/µ dxi ) [(ξi - (∆ξ)xi)Da + dy φLz Re T h )] (xi - qi)(1 - K
(
) Λ(y, xi, xj)
(39)
To apply this model, we consider the palladium membrane WGS microreactor used in fuel microprocessors for micro-fuel-cell applications. Our previous work5,12,20,31,32 has demonstrated the fabrication and complete characterization of a palladium-based membrane microreactor/microseparator used to perform the WGS reaction. In this work, a multicomponent mixture of H2, H2O, CO, CO2, and CH3OH is fed to a microchannel sealed from the top with a palladium membrane to separate H2 from other nonpermeate components. This pure-grade H2 will consequently be supplied to a fuel cell. CO concentration in the permeate stream must not exceed 10 ppm to ensure proper functionality of the proton exchange membrane (PEM) of the fuel cell. To achieve these two objectives, a highly selective palladium membrane is used and a WGS reaction is performed to convert CO into CO2 via the following exothermic reversible reaction in the presence of either Fe2O3/Cr2O3 catalyst at 300 °C or Cu/ZnO/Al2O3 at ∼200 °C:
CO + H2O h CO2 + H2 (∆H298K ) -41.17 kJ/mol) Expressions that describe the reaction kinetics of the WGS reaction have been derived in the literature, based on Langmuir-Hinshelwood kinetics mechanism.33-35 The rate expressions most easy to use were those suggested by Moe,36-38 and Podolski and Kim.39 Amadeo and Laborde’s33 rate expression added an additional term to account for the effect of high pressures at low temperatures. With the catalyst’s density, Fcat, Moe’s expression, which is given in units of mol m-3 min-1, in terms of species concentrations, becomes
9800
Ind. Eng. Chem. Res., Vol. 44, No. 26, 2005
(
r ) 1.85 × 10-5Fcat(RT)2 exp 12.88 1855.5 CCOCH2O[1 - CCO2CH2/(CCOCH2OKe)] (40) T
)
where Ke is the thermodynamic equilibrium constant. Because the inlet stream to the WGS membrane microreactor is assumed to have had reached equilibrium as it emerges from the steam reformer, continuous removal of H2 will shift the reaction equilibrium toward the formation of products (Ke f ∞). In this work, we will use the reaction kinetics suggested by Keiski and co-workers,40 which gave a similar, however, nonelementary, expression for the WGS reaction using Cu/ZnO as a catalyst in terms of the reversibility constant τ (where τ ) [CCO2CH2/(CCOCH2O)]K-1 e ): β (1 - τ) rCO ) k1CRCOCH 2O
(41)
Many researchers assume an Arrhenius-type permeability.16,41-43 However, in light of the permeation mechanism, a different expression for permeability was proposed in our previous work:5,12,30
k hp n h n) (P - P φ s,2
(42a)
kp,0 Ep exp T RT
(42b)
Jp ) k hp )
( )
Our expression dictates an inverse proportionality between the permeability pre-exponential factor and the membrane thickness. This finding is in agreement with what many researchers have observed and is intuitive, from a theoretical standpoint. With thinner membranes, higher fluxes are expected. This expression also shows that permeability does not increase monotonically with temperature. We fabricated our membrane and measured the flux of H2 at different temperatures and upstream pressures, PH2. In view of eq 42 and assuming that the sweep gas always remains fresh (P h ) 0), we fitted our data using the method of least squares to the following expression:
Jp )
[
( )]
Ep kp,0 exp Pn φT RT H2
(43)
From the fit, we determined values of Ep ) 4.0772 kJ/ mol and kp,0 ) 32.5306 mol K m-1 Pa-n s-1. The optimal value of the exponent n, for which the data had a maximum correlation coefficient (R2s ), was n ) 0.489. The latter value is very similar to the value of 0.5 that is customarily used when permeation is believed to be diffusion-limited and, thus, obeys Sievert’s law. 5. Results and Discussion The model that we developed in eq 25 is a nonlinear second-order differential equation initial value problem. An integrator that is forward in space shall be adopted to provide a numerical solution. A Runge-Kutta (RK) integration scheme can be used, because of its proven accuracy and large radius of convergence. To ensure convergence and to avoid numerical solution discontinuities, the spatial step size must lie within the radius of convergence of the RK method. A FORTRAN code was written to perform the calculations. To address the
Figure 2. Mole fractions of the five species and the rate of H2 collected along the length of the membrane microreactor.
accuracy problem, finer meshes were tried, and solution convergence was verified by mesh reduction. The code, which incorporates pressure drop for slipping flows in microdevices,12 can handle nonisothermal operations and uses different mixing rules in numerous subroutines for the evaluation of mixture properties such as viscosity, heat capacities, and species diffusivities. 5.1. Flexibility of the Numerical Solution. Because of the high degree of nonlinearity introduced by incorporating chemical reactions, the general model required a mesh finer than that of the permeation model. Convergence was achieved with 5000 mesh points along the microchannel length of 4 cm. The original second-order model was reduced to a first-order system by utilizing one integral equation for the velocity function. The number of mesh points used to generate a convergent solution for a 4-cm length clearly indicates that the solution would have been substantially more expensive computationally if the model had been a boundary-value problem. Because the primary role of the fuel microprocessor is the production and delivery of H2, a meaningful design would configure the WGS reaction to occur only when the hydrogen concentration falls below the equilibrium limit dictated by the chemical reaction. Starting the reaction before this point will shift the WGS reaction equilibrium toward the consumption of H2 in favor of CO formation, because of the reaction reversibility. This idea suggests using the permeation model until the equilibrium limit is reached, before the general reactive/ separation model is used. The code we developed allows the flexibility of switching between these two different modes of operation. 5.2. Concentration Profiles. Because of velocity slip, a very small drop in the absolute pressure was observed. There was a decrease in the average velocity of the bulk, as a result of mass loss by means of H2 permeation through the membrane. However, fluid density increases with H2 depletion and heavier CO2 formation. Figure 2 shows the cumulative amount of H2 collected, and the concentration profiles of the species along the length of the microchannel. CO is consumed by the WGS reaction along the length of the microchannel, and H2 continues to permeate through the membrane until its partial pressure reaches the low-pressure limit on the other side of the membrane, at which point all H2 has permeated. With larger inlet H2 concentrations and traces of CO, an increase in H2 mole fraction that was due to the WGS was not observed. However, faster reaction rates, lower H2 permeation, or larger inlet CO concentrations lead to H2 buildup in the microchannel. Because of reaction and permeation, the amount of CO2 increases monotoni-
Ind. Eng. Chem. Res., Vol. 44, No. 26, 2005 9801
Figure 3. Spatial variations of numerous dimensionless numbers needed in developing the model relative to their inlet values.
cally. When no reaction is present, all species but the permeate must experience a monotonic increase in their concentration profiles. However, when a chemical reaction besides permeation occurs, only inert and nonpermeating products experience that monotonic increase. Nonpermeating reactants will be consumed by the reaction until the limiting reactant is depleted, after which point their mole fractions will increase, because of permeation. H2O, in this case study, is a nonpermeating reactant, which showed a slight decrease in its mole fraction as the reaction progresses, and then an increase in concentration as permeation continues. 5.3. Dimensionless Numbers and Validation of Assumptions. Figure 3 shows that the decrease in the kinematic viscosity was not comparable to that of the average velocity, hence the decrease in Re. Because isothermal operation was assumed, the concentrationindependent Damko¨hler number (Da), which compares the speed of reaction to that of permeation, remained constant throughout the process. However, Da, which is the modified concentration-dependent Damko¨hler number, which measures the significance of the reaction flux to that of permeation, approaches an infinite value at the end of the microchannel. This suggests the appropriateness of using the no-permeation asymptotic model, as discussed previously in the paper. For isothermal processes, looking at Da alone can be misleading for determining if one of the asymptotic models can be used, instead of the general reactive/separation model. The Knudsen number (Kn), which is the criterion for the slip in microchannel fluid flows, showed a 40% decrease from its inlet value, as shown in Figure 3. This latter result proves that first-order slip models are sufficient to estimate the pressure decrease in this class of devices. The decrease in the Mach number (Ma), as the fluid flows in the microchannel, validates the assumption of flow incompressibility at the inlet of the microchannel, which was assumed during the derivation of the flow model. The Peclet number (Pe), which compares fluid to permeation velocities, remains almost constant until the permeating species almost disappears. The drift velocity of permeation was at least 2 orders of magnitude less than the average axial velocity of the fluid stream in the microchannel. This observation, which is in absolute agreement with the order-ofmagnitude analysis derived from the continuity equation, was anticipated, given our previous work.12 The Prandtl number (Pr) compares viscous momentum effects to the conductive heat transfer, and the Schmidt (Sc) number compares the diffusion of each species in the fluid body. The analysis of the Sc shows that both the viscous- and bulk-diffusion terms are of
Figure 4. Effects of H2 inlet mole fraction on the amount of H2 collected, and the lengths of the membrane and the catalytic reactor. Discrete points are the results from simulation, and the dashed lines are the least-squares fit.
Figure 5. Effect of CO inlet mole fraction on the amount of H2 collected, and the lengths of the membrane and the catalytic reactor. Discrete points are results from simulation, and the dashed lines are the curve fitting using the method of least squares.
the same order of magnitude, thus giving Sc values that are close to unity for all species except H2. Since Sc
