5364
Ind. Eng. Chem. Res. 2005, 44, 5364-5372
Boundary Element Method for Solution of Dispersion Models for Packed Bed Reactors P. A. Ramachandran† Department of Chemical Engineering, Washington University, Campus Box 1198, St. Louis, Missouri 63130
Simulation of heterogeneous models for a packed reactor is a numerically challenging problem due to highly nonlinear coupling of transport, reaction, and the exponential temperature dependence of the rate constants. A system of differential algebraic equations is needed to model the problem. In this paper, the use of the boundary element methods for the numerical solution of the problem is demonstrated. It is shown that integral formulation of the model can be obtained by the elimination of the diffusion-convection operator part of the differential equations by using adjoint formulation. The resulting set of integrals are solved after approximating the dependent variables in each element by a cubic (Hermite) osculation. The efficiency of the method has been demonstrated for a simple exothermic reaction for both pseudo-homogeneous and heterogeneous models, with and without solid-phase conduction. 1. Introduction Mathematical models for reactors such as packed beds are often done using the axial dispersion model.1-3 The models often describe the reactor predictions well and are commonly used as a first-level model to guide the design process. Although the dispersion model has some limitations in the physical description of the process as pointed out in some recent studies,2,4,5 the model is rich in bringing out the potential behavior of the reactor from a practical point of view. Also, it may be noted here that the dynamic behavior of these reactors is not well captured by the dispersion model as pointed out by Parulekar and Ramkrishna,6 but the model offers relatively sound reactor predictions for the steady-state case. The present paper focuses on an efficient numerical procedure for the steady-state solution of the dispersion model. Depending on the relative magnitude of the transport resistances, the packed beds can be modeled by either a pseudo-homogeneous or a heterogeneous model. Detailed comparisons of the pseudo-homogeneous and the heterogeneous models and the range of their validity are presented in a number of studies,3,7,8 and these papers provide the guidelines for the appropriate choice of the type of model. The pseudo-homogeneous model is characterized by a set of second-order boundary valuetype of differential equations, while the heterogeneous model is described by a set of differential-algebraic equations. Due to the interaction of heat and mass transport effects, the solution of the heterogeneous dispersion model has some numerical challenges, especially for highly exothermic reactions. The various techniques for the solution of these classes of problems are discussed now. Finite difference methods are widely used for these problems where the derivative terms are approximated by difference formulas of various orders.9 The orthogonal collocation method developed by Villadsen and Stewart10 is another useful procedure, and this method has been † Tel.: 314-935-6531. Fax: 314-935-7211. E-mail: rama@ wuche.wustl.edu.
applied in many of the earlier studies (for instance, see Young and Finlayson11). Another commonly used solution method is B-spline orthogonal collocation on finite elements,12 and the public domain software programs Colnew and Coldae13 are useful for this purpose. An interesting solution method for second-order differential equations is based on the boundary element method (BEM).14 The method is rooted in the concept of Green’s functions, and the related papers on the use of Green’s function method in chemical reaction engineering15-17 are worth mentioning at this point. The BEM builds on these concepts but is usually applied on subintervals or elements, due to the nonlinear nature of the problem, rather than globally. Compared to other solution methods, the BEM obviates the need to approximate the derivative operators and casts the problem in an integral format. This provides an extra element of a numerical accuracy since the integration is a smoothening process. The method has been applied to diffusion with reaction,18 electrochemical transport,19 and packed bed heat transfer.20 The numerical experiments in these studies show that the number of intervals needed to achieve the same accuracy is considerably less than that for the finite difference method. Further, the nonlinear iterations converge rapidly, even for relatively poor guesses of starting values. The goal of this paper is to illustrate the use of this technique for reactor models with special reference to the dispersion model in a catalytic packed bed. The application of BEM to dispersion model for a single component was shown in an earlier paper by Ramachandran21 and has been extended to multiple systems of differential equations in a recent paper.22 This paper extends the method for heterogeneous dispersion models. The paper also introduces a novel procedure using quasi-linearization for decoupling the particle model from the fluid model. The particle model is not solved separately when using this method and becomes the part of the overall BEM iteration scheme. The extension of the model to include solid-phase dispersion is another novel feature of this paper. It is appropriate at this stage to acknowledge and respect the many contributions of Professor Mike
10.1021/ie049153b CCC: $30.25 © 2005 American Chemical Society Published on Web 03/16/2005
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5365
Dudukovic and wish him many productive years. The pioneering contributions of Mike to reaction engineering are well known. Mike also worked on many novel numerical solution methods, and this paper is an appropriate place to cite some of these papers: Green’s functions in heat transfer;23 Dual series solution methods for diffusion-reaction problems (for instance, Mills and Dudukovic,24); Cumulative concentration concepts for transient problems;25 Moving collocation methods.26 These papers have been well received and have been used by subsequent researchers. The current paper also deals with a novel solution method and is written in the same spirit, and it is one way of honoring Mike’s contributions in this field. The paper is organized in the following manner: Section 2 illustrates the concepts behind the integral formulation and applies this to the pseudo-homogeneous model. Section 3 provides an extension to differential algebraic equations and shows the application to heterogeneous models for packed beds. Section 4 shows the extension for particle phase conduction. For each case, an illustrative example is presented to demonstrate the results of the numerical calculation. Detailed parametric studies related to the reactor modeling and stability are not presented, and no comparison is done against the body of literature in this field since this is outside the main theme of this paper. Some comparison of the numerical results with that predicted by the software Coldae13 is presented to validate the numerical results. Section 5 provides the conclusions and closing remarks. 2. BEM Solution to Dispersion Model 2.1. Problem Statement. The solution methodolgy for the pseudo-homogeneous model is presented in this section, and the extension of to a heterogeneous model is shown in the next section. The prototype problem considered is the same as studied by Ramachandran22 and is as follows: 2 1 d yi dyi ) fi(y1, y2, ‚‚‚, yN) Pei dx2 dx
(1)
Here x represents the dimensionless length (usually 0 to 1). The yi with i ) 1, ‚‚‚, N are the dependent variables. These are the various species concentrations plus the temperature, collectively referred to as species in this paper. The coefficients Pei represent the ratio of characteristic dispersion time to the convection time for the species i. This dimensionless number is, in general, different for different species. The RHS term in eq 1 is a measure of the rate of consumption of the ith species by chemical reaction. For convenience of discussion, eq 1 is multiplied by Pei and rewritten with a coefficient of unity for the second derivative term as shown below, although both representations (eqs 1 or 2) lead to the same result.
d2yi 2
dx
- Pei
dyi ) Pei fi(y1, y2, ‚‚‚, yN) dx
(2)
The boundary conditions are usually of the Danckwerts type, although the limitations of this have been pointed out by Parulekar and Ramkrishna.6 Nevertheless, it is widely used in view of its simplicity and is also used here for illustration of the numerical method. The Danckwerts boundary conditions used here are
stated below: The B.C. at x ) 0 is
[
dyi - Pei yi dx
]
x)0
) -Pei yi0
(3)
where yi0 is the specified value for variable i at the inlet to the reactor (i.e., at x ) 0-). Similarly, B.C. at x ) 1 is represented as
[ ] dyi dx
x)1
)0
(4)
2.2. Adjoint Formulation. The mathematical problem described above can be represented in a compact operator form as
Liyi ) Pei fi(y1, y2, ‚‚‚, yN)
(5)
where the symbol Li is defined as the linear operator for the ith equation and has the following form for the convection-diffusion problem:
Li )
d d2 - Pei dx dx2
(6)
If the operator is of purely diffusion type (e.g., solid conduction operator in Section 4), then
L)
d2 dx2
(7)
Thus eq 5 is a general representation for a linear operator acting on a general forcing function defined by fi. It may be also noted here that the convection-diffusion operator given by eq 6 admits the following asymmetric operator form with a weighting exponential function of the Peclet number.
Li )
d 1 d exp(-Peix) dx exp(-Peix) dx
(
)
(8)
Both the forms are equivalent and the final discretized equations resulting from one formulation is simply a linear combination of the other set. The nonsymmetric form is easier to use from a numerical point of view. Assume that the solution to the eq 5 is needed in the region, x ) a to x ) b (for example, a subdomain or element within 0 to 1) with the boundary values given (or to be calculated) at these end points in terms of either y or the gradient dy/dx or combination thereof. An integral defined as
I)
∫abG(x)L y dx ) 〈G,Ly〉
(9)
is called the inner product of the operator L and an arbitrary function G. Note the symbol 〈*,**〉 for the inner product of two functions as indicated by the last term in eq 9. Taking the inner product on both sides of eq 5 we obtain the adjoint representation which can be written compactly as
〈Gi,L iyi〉 ) 〈Gi,Peifi〉
(10)
which is also the weighted residual representation with Gi being the weighting function.
5366
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
The LHS of eq 10 can be reduced to the following form by integration by parts twice:
〈Gi, L i yi〉 ) Bi + 〈yi, L /i G〉
(11)
which is often referred to as the Green’s formula. The term Bi includes the boundary conditions and is often referred to as the boundary “concomitant” for the ith equation. The specific form of B for the problem studied here will be shown later. The operator L /i in eq 11 is called the adjoint operator for the ith equation. The integral formulations have been done in more purely analytical setting (Arce et al.,16 Ramakrishna and Amundson,28 Dixit and Tavlariades15) and applied to some reactor problems. The procedure applied in these papers is based on the direct use of the Green’s function. Thus the weighting functions are chosen as the Green’s function to the problem defined by the following expression:
L /i G ) -δ(ζ)
(12)
where δ is the Dirac-delta function applied at an interior point ζ in the domain. This leads to an integral equation for the concentration (yi). For an elegant application of this method, the reader should refer to the seminal paper by Arce et al.16 The resulting integral equation was solved by a Picard’s iteration in their work. The procedure adapted here is conceptually similar but somewhat different in details as discussed below. In the BEM formulation presented here, the G functions (also known as weighting functions) are chosen such that
L /i G ) 0
(13)
leading to an integral formulation. Thus homogeneous solutions to the adjoint differential equations are chosen as the weighting functions rather than the actual Green’s function. Also, these solutions are chosen as general solutions to the adjoint equation and do not satisfy any specific boundary conditions, unlike the Green’s function. The overall scheme is more favorable to a direct iterative numerical solution (e.g., Newton iteration), especially for coupled nonlinear equations. The result of using eq 13 in eqs 11 and 10 can be represented in a compact form as
Bi ) 〈Gi,(Pei fi)〉
(14)
which forms the basis for the solution method used in this paper. The specific form of Bi for the convection-diffusion operator can be shown to be the following expression:
[ ] [(
dyi Bi ) G i dx
b
a
-
)]
dGi + PeiGi yi dx
b a
(15)
The adjoint differential equation for this problem for the ith variable is
d 2 Gi dx
2
+ Pei
dGi )0 dx
be noted here that the symmetric formulation given by eq 8 will lead to a symmetric adjoint operator but the inner products will now be defined differently with an exponential weighting function in the integral. Hence, from a numerical point of view both the formulations are equivalent. For the numerical solution, two weighting functions are required since there are four unknowns for each element (concentration and derivative at the two endpoints). Equations for two of these unknowns will come from matching of the neighboring elements or the imposed boundary values at the end points. Hence two additional equations are to be generated for each subinterval for each species. Any two linearly independent solutions to eq 16 can be used in principle as the weighting functions to eliminate the differential operators. However, from a numerical point of view, it is preferable that one of the weighting functions have a value of unity at any given interval. The required weighting functions are chosen as follows:
(16)
Note that the sign change on the convection term indicating that the operator is anti-symmetric. It may
G1i )
1 + exp[-Pei(x - a)] Pei G2i ) 1
(17) (18)
It may also be noted here that the use of the second weighting function (G2i ) 1) is equivalent to imposing an element level mass/heat balance, and hence the current numerical procedure is the mass/heat conserving for each element. The use of these functions in the integral formulation given by eq 14, together with eq 15, provides the following two equations for each subinterval for each species:
-(1/Pei + 1)pia + yia + [1/Pei + exp[-Pei(b - a)]pib - yib )
∫abPeiG1i fi(y1, y2, ‚‚‚, yN) dx
(19)
and
-pia + Peiyia + pib - Peiyib )
∫abPeiG2i fi(y1, y2, ‚‚‚, yN) dx
(20)
where i is the index for each differential equation. Note that G2i ) 1 in eq 20 but is shown there for generalization. The variable pia is defined here as gradient of yi (i.e., dyi/dx) at x ) a, and pbi is the gradient x ) b. Also yia and yib are the values of yi at the positions a and b respectively. The set of eqs 19 and 20 represent the integral representation to the differential equation corresponding to the variable i. If the total number of elements or subintervals used are M, then there are M + 1 nodes and hence 2N*(M + 1) total unknowns, viz., yi and pi for each species at each node. Equations 19 and 20 are applied for all these M elements to generate 2N*M equations. These equations together with 2N boundary conditions provide all the necessary equations for the 2N*(M + 1) variables. The final assembled set of equations have a banded structure and hence the method can handle a fairly large number of reaction species. 2.3. Numerical Representation. The integral representation given by eqs 19 and 20 are the exact representation of the problem. In particular, eq 20 is a
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5367
statement of the element level mass/heat balance, and hence the solution is also automatically mass conserving. However, some numerical approximations are needed to evaluate the integrals on the RHS of these equations and to develop an iterative procedure to arrive at a converged solution. To evaluate the integrals, the profile of yi for any particular element is approximated by an osculating polynomial, as done in earlier papers.14,19,20 Thus local or element level representation for yi is
yi ) φ1(b - a)pia + φ2yia + φ3(b - a)pib + φ4yib (21) where the osculating (Hermite) polynomials φi (i ) 1 to 4) are defined as follows:
Table 1. Coefficients Needed in Equation 31a hki ) ∫baGki dx Rk1i ) (b - a)∫baGkiφ1(η) dx Rk2i ) ∫baGkiφ2(η) dx Rk3i ) (b - a)∫baGkiφ3(η) dx Rk4i ) ∫baGkiφ4(η) dx a k ) 1 and 2 corresponding to the two weighting functions. Also, the coefficients need to be computed for each i since the weighting functions Gki are different for each i.
Table 2. Solution Values for the Base Case distance
concentration
temperature
0.0 0.2 0.4 0.6 0.8 1.0
0.9982 0.8765 0.6946 0.4014 8.43E-02 7.61E-03
1.03E-03 2.56E-02 6.26E-02 0.1221 0.1842 0.1984
φ1 ) η - 2η2 + η3
(22)
φ2 ) 1 - 3η2 + 2η3
(23)
the above expression leads to
φ3 ) -η2 + η3
(24)
Iki ) Peikoi
φ4 ) 3η2 - 2η3
(25)
∫abGki dx + N
Pei
∫abGki( ∑ k1im[φ1(b - a)pma + φ2yma + m)1
with η defined as the dimensionless distance variable within each subinterval: η ) (x - a)/(b - a). Note that the profile given by eq 21 not only fits the values at the end points of the subinterval but also has the correct slope at these points. For the purpose of the iterative solution, the rate term can be adequately represented by a local linearized version. This scheme is essentially a variation of the Newton-Raphson iteration. Using linear terms of the Taylor series, the nonlinear rate fi can be represented as N
fi ≈ koi +
∑ k1imym
(26)
m)1
where koi is the constant part of the Taylor expansion of fi defined as N
koi ) f /i -
( ) dfi
y/m ∑ dy m)1
*
(27)
m
and k1im is the linear term and corresponds to the Jacobian matrix of the rate function
( )
dfi * k1im ) dym
(28)
where the * denotes the elemental average values at the current iteration. An average value of ym over the subinterval can be used to evaluate these linearized coefficients in eqs 27 and 28. Thus y/m ) (yma + ymb)/2. Also note that these coefficients change at every iteration. The integral term in eqs 19 and 20 (represented here as Iki; k ) 1 and 2) can then be expressed as N
∫abGki(koi + ∑ k1imym) dx
Iki ) Pei
(29)
m)1
Substitution of the osculation approximation for ym in
φ3(b - a)pmb + φ4ymb]) dx (30) The above equation can be written in the following form by defining a set of element level integrals: M
Iki ) Peikoihki + Pei
∑ (k1im[Rk1ipma + Rk2iyma + m)1 Rk3ipmb + Rk4iymb]) (31)
Equation 31 provides the numerical representation for the integral terms in eqs 19 and 20. The coefficients in the above equation are the inner products of the weighting functions and the osculation polynomials and are defined in Table 1. These coefficients do not change with the iterations and can be calculated (for a given set of subintervals) at the beginning of the numerical procedure. 2.4. Test Example. A system of two differential equations representing a single exothermic reaction was chosen as the first test example. The variable y1 here is the concentration and y2 is the dimensionless excess temperature and the rate functions are
f1 ) Da y1 exp
( ) ( ) γy2 1 + y2
f2 ) -BDa y1 exp
γy2 1 + y2
(32)
(33)
where the dimensionless parameters are standard. Da is the reaction Damkohler number (ratio of convection time to reaction time), γ is the dimensionless Arrhenius energy, and B is the dimensionless heat of reaction. The parameters chosen were as follows: Pe1 ) 300; Pe2 ) 100; Da ) 0.5; B ) 0.2 (base case); and γ ) 20. The results for this case generated with equally spaced 10 subintervals are shown in Table 2 and agree well with those obtained by a finer set of 41 subintervals. Results were also validated using the program Colnew. The parametric effect of heat of reaction is then studied in order to validate the numerical methods for
5368
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
Table 3. Solution Values for B ) 0.3 with 10 Intervals; Other Parameters Same as the Base Case distance
concentration
temperature
0.0 0.2 0.4 0.6 0.8 1.0
0.9983 0.8533 0.4051 1.0E-06 1.0E-06 1.0E-06
1.58E-03 4.61E-02 0.1908 0.3005 0.3000 0.2999
a problem with steep gradients. The B parameter was changed to 0.3 while all the other parameters were kept the same as base case. Again, 10 subintervals were used initially and the results are shown in Table 3. Note that the negative concentrations in the iterations were set as 1 × 10-6 in this simulation. Small negative concentrations and slight wiggles in temperature are observed, and these disappear by adaptive remeshing shown later. The direct solution with Colnew did not converge for this case since it is sensitive to the starting values. It converged only after a continuation starting at a lower value of the B parameter and progressively increasing this parameter to 0.3. In contrast, the BEM method appears to have a wider window of starting values for which the convergence is obtained. The temperature gradients are also of importance in stability analysis and as a guide to optimum location of mesh intervals. The gradients are directly obtained as a part of the solution in the boundary element method. The results for the gradients as a function of axial distance are shown in Figure 1 for B ) 0.3 and the sharp peak at x ) 0.4 should be noted. The second derivatives calculated from the first derivative information by finite differences is shown in Figure 2. Note that the second derivative changes sharply over a short distance going from a large positive value to a large negative value. The finite difference (FD) approximations have difficulties in accurately capturing this behavior. This may be the reason for the poorer convergence of these methods and the need for very fine meshes in FD near the ignition point (in this case near 0.4). The present method overcomes these problems. The collocation method (e.g., Colnew) will also require a large number of subdomains in this region. The use of the gradient information as a tool for adaptive improvement is now illustrated. The result of Figure 1 suggests that the mesh points should be relocated near x ) 0.4 (for the case of B ) 0.3). A simple way of generating nonuniformly spaced intervals is to use the roots of Jacobi polynomials as the mesh points. These polynomials were introduced into the chemical engineering literature by Villadsen and Michelson.29 They are orthogonal with respect to the weighting function of xR(1 - x)β. If R and β are both zero, then roots are symmetric and correspond to the standard Gauss points. If R is chosen larger than β, the roots are more closely placed near x ) 0 than x ) 1 and vice versa. The values of the first derivatives can be used to find the location where the nodes should be more dense. For the above example (for B ) 0.3) we find the nodes should be near x ) 0.4 which is close to the ignition point for the reaction. The Jacobi polynomial with R ) 0 and β ) 4 is chosen for the region 0 < x < 0.4; similarly, R ) 4 and β ) 0 were chosen for 0.4 < x < 1. This automatically provides an adaptive mesh near the ignition point. Figure 3. shows the concentration and the temperature profiles generated by this nonuniform mesh with only
Figure 1. First derivative of temperature for homogeneous model B ) 0.3. Other parameters same as base case.
Figure 2. Second derivative of temperature for B ) 0.3 showing widely differing values in the second derivative near the reaction front.
Figure 3. Pseudo-homogeneous model: Concentration and temperature profiles for B ) 0.3 with uneven subintervals.
14 subintervals. No wiggles or negative concentrations were obtained as shown in Figure 3. 3. Heterogeneous Model 3.1. Model Equations. For the heterogeneous model, the fluid phase mass balances and the heat equation
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5369
have the following general form:
d2yi 2
dx
- Pei
dyi ) PeiRi(yi - zi) dx
Table 4. Solution Values for Heterogeneous Model for the Base Case; Direct Solution with the Starting Values of Zero
(34)
where Ri represents the dimensionless mass/heat transfer coefficient for species i. The variables zi are the corresponding particle variables, viz., particle surface concentrations and particle temperature. It is convenient to write the above equation in a vector-matrix form as
dy b y d 2b -P ˜e ) P ˜ eR˜ (y b-b z) dx dx2
(35)
where R˜ is a diagonal matrix with the values of Ri on the ith diagonal and zero elsewhere. Similarly, P ˜e is a diagonal matrix with the values of Pei on the diagonals. The particle models are generally of the following form:
Ri(yi - zi) ) Fi(z1, z2, ‚‚‚)
(36)
where Fi is the rate of consumption of the species i within the catalyst particle (based on unit reactor volume). The equation can be represented in vector matrix form as
R˜ (y b-b z) ) F B
(37)
The system is thus a differential-algebraic set given by eqs 35 and 37. 3.2. Solution Procedure. The RHS functions of eq 37 is linearized at the current level of iteration and is represented as
˜b z F B)K B0 + K
(38)
where K B 0 is the constant term in the linearization and K ˜ is the Jacobian matrix of the rate functions. These coefficients are defined below and are evaluated using values at the current level of iteration (denoted by superscript *) for the solid phase N
K0i ) F/i -
( ) dFi
z/m ∑ dz m)1
*
(39)
m
and Kim is the linear term and corresponds to the Jacobian matrix of the rate function
Kim )
( )
dFi * dzm
(40)
Substitution of 38 in eq 37 yields
˜b z R˜ (y b-b z) ) K B0 + K
(41)
which can be solved for the z vector
b z ) (K ˜ + R˜ )-1R˜ b y - (K ˜ + R˜ )-1K B0
(42)
and is the partial solution (one step of Newton’s iteration) for the particle model. In the current solution scheme this solution is directly used in the reactor
x
y1
z1
y2
z2
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.9978 0.9217 0.7627 0.1819 0.0730 0.0293 0.0117 0.0047 0.0019 0.0008 0.0003
0.9815 0.8998 0.7009 0.1395 0.0559 0.0224 0.0090 0.0036 0.0014 0.0006 0.0002
0.0014 0.0169 0.0588 0.1559 0.1765 0.1848 0.1881 0.1894 0.1900 0.1902 0.1903
0.0144 0.0344 0.1083 0.1899 0.1902 0.1903 0.1903 0.1903 0.1903 0.1903 0.1903
model. Thus the Bf vector needed in the fluid phase balance
Bf ) R˜ (y b-b z)
(43)
is obtained directly in a linear form as
Bf ) B k 0 + k˜ 1b y
(44)
B k 0 ) R˜ (K ˜ + R˜ )-1K B0
(45)
˜ + R˜ )-1R˜ ] k˜ 1 ) R˜ [I˜ - (K
(46)
where
and
where I˜ is the identity matrix. The fluid model then reduces to a linear form containing only the vector b y and can be readily solved by the procedure given in Section 2. The solution procedure used can then be summarized as follows. (i) At the current level of iteration for the gas phase, the values of the concentrations in the solid phase (and temperature) are used to obtain the linearized ‘rate’ terms for the solid. (ii) These are then used to compute the coefficients needed in the gas phase balance by eqs 45 and 46. (iii) The concentrations in the solid are also updated at this point using eq 42. (iv) The BEM formulation for gas phase (which is linear now) is now solved to get updated values for the gas phase variables. The procedure is repeated until convergence is attained, which usually takes six to seven iterations. 3.3. Test Example. A single exothermic reaction was again chosen as the test example. The particle functions F1 and F2 were chosen the same as were used for the pseudo-homogeneous model test example 2.3 given by eqs 32 and 33. The same parameters were used. The additional parameters needed are the mass and heat transfer coefficients, and these were set as R1 ) 40 for mass transfer and R2 ) 10 for heat transfer. Results for this set of parameters show a multiple steady-state behavior that is not observed for the pseudo-steady-state model. The results are shown in Table 4 for one of these steady states. These results are obtained by direct iteration with a starting value of zero temperature. The solution for the same set of parameters was also obtained by starting with an R2 of 5 and continuing the calculations by progressively increasing this parameter to 10. This procedure yielded the upper ignited steady state, and the resulting solution is shown in Table 5. The particle temperature at x ) 0 is much higher in the ignited state and even larger than the adiabatic temperature rise of 0.2 for this case.
5370
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
Table 5. Solution Values for the Heterogeneous Model for the Same Case as in Table 4 Obtained by Continuation from a Smaller Value of the Heat Transfer Parameter x
y1
z1
y2
z2
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.9210 0.2430 0.1114 0.0495 0.0218 0.0095 0.0042 0.0018 0.0008 0.0003 0.0002
0.0126 0.1969 0.0886 0.0392 0.0172 0.0075 0.0033 0.0014 0.0006 0.0003 0.0001
0.0368 0.1348 0.1599 0.1717 0.1770 0.1793 0.1803 0.1807 0.1809 0.1810 0.1811
0.7635 0.1716 0.1781 0.1800 0.1806 0.1809 0.1810 0.1811 0.1811 0.1811 0.1811
The profile of the particle temperature at x ) 0 as a function of the parameter R2 was also generated, and the results are shown in Figure 4. The lower curve (marker *) was obtained by starting from R2 ) 20 and using the converged solution for each case as the starting value and progressively decreasing the value of R2. This yields the lower steady state with the particle temperature at x ) 0 close to the inlet gas temperature. The upper curve in Figure 4 (marker +) shows the results of a continuation with a progressive increase in the heat transport parameter. The staring value was R2 of 5. The results obtained show an ignited behavior for certain ranges of R2, and in this range considerable temperature difference between the gas and solid is observed for this case. Detailed continuation studies and bifurcation diagrams are not studied here since the purpose of the paper is to demonstrate the solution method and also to indicate that the method is capable of providing multiple solutions in certain range of parameters. The combination of BEM with continuation methods such as DERPER30 is useful for such problems, and the implementation will be reported in a later paper.31 It may be noted that the orthogonal collocation method with the bifurcation software AUTO has been used to do a bifurcation analysis of the dispersion model in a recent study by Liu and Jacobson.32
4. Particle Conduction Model The direct conduction between the solid phase (at the contact points) may be an important factor, especially for the case of large solid-phase temperature gradients and for high-temperature reactions. The model equation for the solid phase is now modified to include the conduction term containing the second derivative of the solid temperature. The BEM solution method to this problem is discussed here. The problem considered here is represented as follows. The gas-phase mass balance equation is the same as eq 34 and is reproduced here for convenience.
d2y1 dx
2
- Pe1
dy1 ) Pe1R1(y1 - z1) dx
(47)
The gas phase heat balance has also the similar form and is
d2y2 2
dx
- Pe2
dy2 ) Pe2R2(y2 - y3) dx
(48)
Here we use the notation y3 to represent the particle temperature, and the differential equation for this is
d2y3 2
dx
) Pe3R2(y3 - y2) - Pe3BDa z1 exp
( ) γy3 1 + y3
(49)
The Danckwerts boundary conditions are used for y1 and y2, while for y3 we see the no flux boundary condition at both x ) 0 and x ) 1. Finally, the mass transport to the particle is given by the algebraic equation
R1(y1 - z1) ) Da z1 exp
( ) γy3 1 + y3
(50)
The system now consists of three differential equations and one algebraic equation. The algebraic equation in general needs an iteative solution but for this case (first order reaction) the equation is simple and directly solved as
z1 )
y1
( )
1 + (Da/R1) exp
γy3 1 + y3
(51)
Thus we have three differential equations with mixed operators. Two equations for gas phase have the convection-diffusion operator while the solid phase equation has now only the pure diffusion operator. The weighting function for the solid phase heat equation is chosen differently now. The adjoint equation for the diffusion operator is also the diffusion operator since this is a symmetric operator. The G functions are therefore linear combinations of 1 and x. Convenient choices for eq 49 are therefore
Figure 4. Heterogeneous model: Plot of particle temperature at the reactor inlet as a function of the heat transfer parameter R2. * is the curve obtained by doing a decreasing continuation on the R2 parameter. + is the curve doing an increasing continuation. The plot shows that the particle temperature can have multiple values for a range of R2.
G13 ) x - a
(52)
G23 ) 1
(53)
and
The boundary term, B (see eq 14), needed in the integral
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005 5371 Table 6. Solution Values for the Test Example 3.3 Including the Effect of Solid Conduction distance
gas conc.
gas temp.
particle temp.
0.0 0.2 0.4 0.6 0.8 1.0
0.927 0.049 0.006 0.001 0.000 0.000
0.046 0.190 0.199 0.200 0.200 0.200
0.327 0.204 0.200 0.200 0.200 0.200
formulation for the variable 3 is now given as
[
] [ ]
dy3 B3 ) G3(x) dx
b
a
-
dG3 y dx 3
b a
(54)
With these modifications the solution procedure is similar to that in section 2.3. An illustrative result for test example 3.3 is shown in Table 6. Here we use the same parameters as the base case in section 3.3, except for the heat Peclet number. These are now split between the gas and solid as Pe2 ) 50 and Pe3 ) 50. The results in Table 6 show that the solid temperatures are smaller than the values in the model neglecting the effect of solid conduction. This indicates the importance of including the solid phase conduction, particularly for thermal runaway and stability analysis. This also indicates the need for more descriptive mesoscale models for packed beds to generate the parameters Pe2 and Pe3 in the dispersion model. No multiple solutions were found for this case, but a detailed bifurcation analysis would be needed to confirm this. More rigorous models for particle conduction were discussed by Balakotaiah and Dommeti.8 These involve even higher derivatives of temperature and hence are computationally difficult by traditional finite difference methods. The method proposed here eliminates the need for approximations for the derivatives and is hence suitable for such complex problems. These refinements are not studied here and could form the basis for future work. 5. Concluding Remarks Boundary integral methods have been shown to be useful for the solution of the dispersion model. The method uses the adjoint representation and hence eliminates the need to approximate the first and second order derivatives. This leads to an increased accuracy with relatively fewer number of elements. Also, the solution converges rapidly even for starting values that are not quite close to the exact solution. The method has built-in element level mass/heat conservation properties, which also lead to an improved accuracy. The first derivatives at the various points are also directly obtained as a part of the solution. This information is useful to assess the runaway characteristics and parametric sensitivity studies. The first derivative information is also useful in an adaptive element placement, which can be done easily by using the roots of the Jacobi polynomials. For the heterogeneous model, an iterative scheme is developed based on quasi-linearization of the particle model by which the particle model is not separately solved but becomes a part of the overall iterative scheme. One test example examined indicated multiple steady states.
The method was also extended to the heterogeneous model with solid conduction effects. This is an interesting case since the operators for the various equations are different. Thus, weighting functions for each equation are different and depend on the linear operator for that equation. This is in contrast to the finite element method where the weighting functions are all same for all the equations and correspond to the interpolating functions. The test case for solid conduction shows that the solid temperatures are somewhat smeared with the solid conduction and the relatively large temperature at the ignited state is reduced. This shows the importance of including solid conduction for highly exothermic reactions, especially for thermal stability and control studies. The method was developed for the constant velocity case and constant values of the transfer coefficient but can be modified to include the varying effects as well. Since each element is separate and coupled to the adjacent element only through the continuity of concentration and boundary fluxes, the various transport parameters can be different for each element. Thus variable properties can easily be included in the model although we have not shown any test example for these cases. Each element acts as a small packed bed reactor. Similarly, variable velocity in the convection term can be accommodated as a perturbation term and will contribute an extra integral term, and the procedure can be extended for this case as well. Literature Cited (1) Rase, H. F. Fixed Bed Reactor Design and Diagnostics; Butterworths: Boston, 1990. (2) Iordanidis. Mathematical Modeling of Catalytic Fixed Bed Reactors. Doctoral thesis, Twente University, 2002. (3) Dommeti, S, M.; Balakotaiah, V.; West, D. H. Analytical Criteria for Validity of Pseudohomogeneous Models of Packed Bed Catalytic Reactors. Ind. Eng. Chem. Res. 1999, 38, 767-777. (4) Kronberg, A. E.; Westerterp, K. R, Nonequilibrium effects in fixed bed interstitial fluid dispersion. Chem. Eng. Sci. 1999, 54, 3977-3993. (5) Chakraborthy, S.; Balakotaiah, V. Low-dimensional models for describing mixing effects in laminar flow tubular reactors. Chem. Eng. Sci. 2002, 57, 2545-2564. (6) Parulekar, S.; Ramkrishna, D. Analysis of axially dispersed systems with general boundary conditions. Chem. Eng. Sci. 1984, 39, 1571-1579. (7) Ramkrishna D.; Arce, P. Can Pseudo-homogenoeus reactor models be valid? Chem. Eng. Sci. 1989, 44, 1949-1969. (8) Balakotaiah, V.; Dommeti, S, M. Effective models for packed bed catalytic reactors. Chem. Eng. Sci. 1999, 54, 1621-1628. (9) Schiesser, W. E. Computational Mathematics in Engineering and Applied Science: ODEs, DAEs and PDEs; CRC Press: Boca Raton, 1994. (10) Villadsen, J. V.; Stewart W. E. Orthogonal Collocation Methods in Chemical Engineering. Chem. Eng. Sci. 1967, 22, 1483-1496. (11) Young, L. C.; Finalyson, B. A. Axial dispersion in nonisothermal packed bed chemical reactors. Ind. Eng. Chem. Fundam. 1973, 12, 612-618. (12) Ascher, U. J.; Christiansen, J.; Russell, R. D. COLSYS: A collocation software for boundary value problems, Codes for Boundary value problems. Chides, et al., Eds.; Lecture notes in Computer Science 76; Springer-Verlag: Berlin, 1978. (13) Dension, K. S.; Hamrin, C. E.; Fairweather, G. Solution of Boundary Value Problems Using the software Package DS04AD and COLSYS, Chem. Eng. Commun. 1990, 90, 103-110. (14) Ramachandran, P. A. Boundary Element Methods in Transport Phenomena; Computational Mechanics Publications: Southampton, U.K., 1993. (15) Dixit, R. S.; Tavlarides, L. L. Integral method of analysis of Fischer-Tropsch synthesis reaction in a catalyst pellet. Chem. Eng. Sci. 1982, 37, 530-544.
5372
Ind. Eng. Chem. Res., Vol. 44, No. 14, 2005
(16) Arce, P. A.; Cassano, A.; Irazoqui, H. Tubular reactor with laminar flow regime: An integral equation approach. Comput. Chem. Eng. 1988, 12, 1103-1114. (17) Ramkrishna, D.; Arce, P. Self-adjoint operators of transport in interacting solid-fluid systems. Chem. Eng. Sci. 1988, 43, 933944. (18) Ramachandran, P. A. Application of the boundary element method to nonlinear diffusion with reaction problems. Int. J. Numerical Methods in Eng. 1990, 29, 1021-1031. (19) Barbero, A. J.; Mafe, S.; Ramirez, P. Application of the boundary element method to convective electrodiffusion problems in charged membranes. Electrochim. Acta 1994, 39, 20312035. (20) Lesage, F.; Latifi, M. A.; Midoux, N. Boundary element method in modeling of momentum transport and liquid-to-wall mass transfer in a packed bed reactor. Chem. Eng. Sci. 2000, 55, 455-460. (21) Ramachandran, P. A. Boundary integral solution method for axial dispersion with reaction problem. Chem. Eng. J. 1990, 45, 46-54. (22) Ramachandran, P. A. Boundary integral solution method for multi-species 1-D convection-diffusion-reaction problems. Comput. Math. Appl., in press. (23) Dudukovic, M. P. Heating of an insulated slab by nonuniform heat generation. Chem. Eng. Sci. 1981, 5, 109-110. (24) Mills, P. L.; Dudukovic, M. P. Application of the method of weighted residuals to mixed boundary problems: Dual series. Chem. Eng. Sci. 1980, 35, 1557-1570. (25) Dudukovic, M. P.; Lamba, H. S. Solution of moving boundary problems by orthogonal collocation. Chem. Eng. Sci. 1978, 33, 303-309.
(26) Ramachandran, P. A.; Dudukovic, M. P. Moving finite element collocation method for transient problems with steep gradients. Chem. Eng. Sci. 1984, 35, 1321-1324. (27) Arce, P. A.; Locke, B. R.; Trigatti, I. M. B. An integralspectral approach for reacting Poiseuille Flows. A. I. Ch. E. J. 1996, 42, 23-41. (28) Ramkrishna, D.; Amundson, N. R. Linear Operator Methods in Chemical Engineering; Prentice-Hall: Englewood Cliffs, New Jersey, 1985. (29) Villadsen. J.; Michelson, M. L. Solution of Differential Equation Models by Polynomial Approximation; Prentice-Hall: Englewood Cliffs, New Jersey, 1978. (30) Kubicek, M.; Marek, M. Computational Methods in Bifurcation Theory and Dissipative Structures; Springer-Verlag: New York, 1983. (31) Ramaswamy, R. C.; Ramachandran, P. A. Multiple steadystate calculations for diffusion problems using BEM and DERPER. Manuscript in preparation. (32) Liu, Y.; Jacobson, E. W. On the use of reduced order models in bifurcation analysis of distributed parameter systems. Chem. Eng. Sci. 2004, 28, 161-169.
Received for review September 5, 2004 Revised manuscript received December 15, 2004 Accepted December 17, 2004 IE049153B
