A modified orthogonal collocation method for reaction diffusion problems

PROCESS SYSTEMS ENGINEERING

A modified orthogonal collocation method for reaction diffusion problems

M. Soliman^{I, II}; Y. Al-Zeghayer^I; A. AjbarIII, ^*

^IDepartment of Chemical Engineering, College of Engineering, King Saudi University, P.O. Box 800, Riyadh, 11421, Saudi Arabia. E-mail: [e-mail]

^IIDepartment of Chemical Engineering, The British University in Egypt, El-Shorouk City, Egypt. E-mail: [e-mail]

^IIIDepartment of Chemical Engineering, King Saud University, P.O. Box 800 Riyadh 11421, Saudi Arabia. Phone: ++ 966-1-4676843, Fax: ++ 966-1478770. E-mail: [e-mail]

ABSTRACT

A low-order collocation method is often useful in revealing the main features such as concentration and temperature profiles and the effectiveness factor for porous catalyst particles. Two modifications are introduced in this paper to make the method more efficient. The first modification is to add an extra collocation point at the center of the particle. It is shown that such extra point introduces a single variable non-linear equation to be solved after obtaining the standard collocation method solution. In the second modification, the polynomial solution obtained from the application of the orthogonal collocation method is transformed to a rational function form. These two modifications are applied to specific examples and it is shown that they can improve the performance of collocation methods in general and the one-point collocation method in particular.

Keywords: Orthogonal collocation; Jacobi polynomials; Interpolation; Reaction diffusion; Catalyst particle.

INTRODUCTION

In the method of weighted residuals, we seek a solution for a differential equation in terms of a polynomial with unknown coefficients to be determined such that certain criteria are satisfied. These criteria are usually chosen such that the integrals of the weighted residual, which is obtained by substituting the assumed solution in the differential equation written with one of its sides as zero, are zeros. Depending on the weight, we obtain methods like the Galerkin method, least square methods, the method of moments, and collocation methods. In collocation methods, we would like to have the residual to be zeros at particular points called collocation points. If the collocation points are the zeros of orthogonal polynomials, the method is called orthogonal collocation. In all methods of weighted residuals, a non-linear ordinary differential equation is approximated by a set of non-linear algebraic equations.

The orthogonal collocation method was developed by Villadsen and Stewart (1967). It was then subsequently studied and applied to many chemical engineering problems by several authors (Finalyson, 1972; Villadsen and Michelsen, 1978; Finalyson, 1980; Torres et al., 2000; Biscaia Junior et al., 2001).

Recently, rational methods were also developed (Berrut et al., 2005). They are based on approximating the solution for a differential equation by a rational function, which is usually more accurate than a straightforward polynomial. Rather than starting with a rational function whose parameters are to be determined, we use a different approach in this paper. Having obtained a polynomial solution, we convert this solution into a rational function.

Solution of catalytic reactor models requires the solution of a set of differential equations representing the phenomena of diffusion with chemical reaction inside the pores of the catalyst particles at each location along the reactor. It is clear that such calculations would require excessive computation time. It is our objective in this paper to make the method of orthogonal collocation more efficient in order to reduce the computational load.

In the next sections we introduce Jacobi polynomials, and then we introduce the standard method of orthogonal collocation. This is followed by presenting the modifications. Then we derive formulae for the effectiveness factor calculation. Finally we apply the method to some problems of diffusion and reaction in catalyst particles.

JACOBI POLYNOMIALS

Usually orthogonal polynomials are defined on the closed interval [-1, 1]. However by simple linear transformation, we can change the definition to any other closed interval. For the purpose of this study we would like to define the polynomials in the interval [0, 1]. In this case they are called "shifted" polynomials. However, we will drop the word "shifted" for brevity.

The Jacobi polynomials of degree n are defined such that they satisfy the orthogonality conditions (Villadsen and Stewart, 1967)

where w(x) is the weighting function for the orthogonality conditions and Cn is a constant. For Jacobi polynomials,

Thus α, β are the indices of the weight function w(x). We will drop the superscripts α, β for brevity.

Jacobi polynomials satisfy the differential equation:

For the case of α = 1 which is of interest in this work, this equation can be written as,

The first member of the Jacobi polynomials, P₀, is of course a constant; thus from Equation (1) for α = 1, we have:

MATHEMATICAL FORMULATION AND SOLUTION BY THE STANDARD COLLOCATION METHOD

Mathematical Formulation

Consider a reaction with a dimensionless rate R(u), where u is the dimensionless concentration. The reaction takes place isothermally inside a catalyst particle with no external resistance to mass transfer. We may write the describing equation as:

with the boundary conditions:

Here u is the dimensionless concentration of the key component in the reaction, r is the dimensionless space variable in the catalyst, is the Thiele modulus and s is the shape factor of the catalyst (s =0 for an infinite slab, s =1 for an infinite cylinder and s = 2 for a sphere). The dimensionless reaction rate R(u) is normalized with respect to the reaction rate at the external surface of the catalyst and, hence, at the external surface r=1, u = 1 and R (1) = 1. The effectiveness factor η is given by:

Since the boundary condition (9) implies that the solution will be symmetric with respect to (r), we introduce the transformation:

so that the problem is changed to:

and

The effectiveness factor is given by:

Except for a few cases, such as isothermal first-order and zero-order reactions, the analytical solution of the boundary value problem, Eqs. (7)-(9) or Eqs. (12), (13) is, in general, not feasible, and the problem can only be solved numerically. In the case when the effectiveness factor is calculated repeatedly, such as in the simulation of packed bed catalytic reactors, we need a fast and efficient numerical method. A simple, convenient and easy approach to solve such problems is through the use of the orthogonal collocation method. It is the purpose of this paper to introduce two modifications to the collocation method which may yield somewhat higher accuracy for the concentration profile and the effectiveness factor than the standard collocation method.

Collocation Solution Representation

In the orthogonal collocation method, we represent the solution as a polynomial with unknown coefficients to be determined such that the differential equation is satisfied at certain points. A better way of writing a polynomial approximation of a function in terms of its values (ordinates) at certain points is through the Lagrange interpolation formula. It has many advantages compared to the straightforward polynomial. First the solution is obtained directly at collocation points. Secondly, when we have a non-linear differential equation, we could have a very good initial guess for the ordinates ui's which lie between 0 and 1, whereas we do not know the range of the parameters in the straightforward polynomial. Thirdly and most important, we can easily derive formulae for the derivatives using the Lagrange interpolation formula. Fourthly, the same Lagrange formula can then be used to obtain the value of u at any point x [0, 1].

The Lagrange interpolation formula is given by:

where

are the Lagrange polynomials.

From Equation (15), the first and second order derivatives of any function u(x) expressed in terms of Lagrange polynomials and at a collocation point x = xi can be written as:

where

An efficient method for calculating the elements a_ij and b_ij is given in the book of Villadsen and Michelsen (1978).

We can also present quadrature Radau formulae (Villadsen and Stewart (1967)) for the effectiveness factor in the form:

The weights of the quadrature are obtained by the procedure described by Villadsen and Michelsen (1978). This formula uses one extra point in x=1, and the zeros of the Lagrange polynomial are calculated using (α+1) instead of only α (Villadsen and Michelsen, 1978).

The collocation points are the zeros of the Jacobi polynomials suitable for the geometry of the catalyst particle. In our case, we use the Jacobi polynomials indices α = 1, β = (s-1)/2

Application of One-Point Standard Collocation to a First-Order Reaction

For a first-order polynomial with the function value given as u(1) at x = 1 and u (x₁) at x = x₁, the Lagrange interpolation formula is given by:

This formula satisfies the values of u at x = 1, and at x = x₁.

Now if we choose x₁ to be the zero of a Jacobi orthogonal polynomial P₁(x) such that α = 1, β = (s-1)/2, we obtain:

For a first-order reaction, Equation (12) becomes:

and u(x) becomes:

One can show that this solution is accurate up to terms containing at x = x₁ and up to elsewhere. In addition, the effectiveness factor, given by the following Equation (26), is accurate up to

MODIFIED COLLOCATION METHOD

Two modifications to the standard collocation method are now suggested to improve its performance.

First Modification

One can show that, if we collocate the equations at an extra point in addition to those obtained as the zeros of the proper orthogonal polynomial, we will have the ordinate of this extra point appearing only in the collocation equation of this point. This ordinate will not appear or affect the other equations.

Let us add an extra collocation x₀ point, preferably at x = 0, and define u_e such that

We notice that the term multiplying is given by:

which, according to Equation (5), when it is subjected to the diffusion differential operator (Equation (12)) will give a term that contains a Jacobi polynomial whose value is zero at the collocation points. Thus, the collocation equations and their solution do not change with the addition of the extra point. This extra point will also not affect the ordinates at the collocation points. Now the collocation at this extra point will give a non-linear equation in a single unknown variable which can be solved after solving the standard collocation equations. This extra point will also not affect the effectiveness factor.

Application of One Point Collocation First Modification (mod 1) to a First-Order Reaction

Let us write the Lagrange interpolation formula in terms of ordinates at one interior collocation point x₁ besides at the boundary x = 1. We have:

Let us define another solution u_e such that its Lagrange interpolation formula in terms of ordinates at an interior collocation point x₁, besides at the boundary x = 1 and any point x₀, is:

Now if x = x₁, we note that and both

at x = x₁ and the integral do not contain a term of u_e (x₀). This means that, whatever the value of u_e (x₀), u(x₁) and the effectiveness factor will not change. This is the major strength of the collocation method as applied to the catalyst particle problem. The choice of x₀, and thus u_e (x₀), will only affect the profile of u_e (x). Now for an arbitrary x₀, we could write from Equation (12):

Thus, after calculating u(x₁) from Equation (25), we calculate u_e (x₀) from Equation (31) to obtain:

and

Now x₀ can be chosen to improve the accuracy at a certain point. x₀ can be chosen as zero so that we obtain a higher accuracy for the value of (u) at the center of the catalyst particle. In Soliman (1988) the choice was made.

In summary, the modified one point collocation requires the solution of two algebraic equations sequentially and it makes the whole profile exact for terms up to , whereas in the original one point collocation the solution is exact for terms up to only at the collocation point and up to elsewhere.

Second Modification

Let us form a rational function from the polynomial solution , and such that:

where

and C is a constant to be estimated such that:

Thus

This choice makes u_r(x) equal to u_e(x) at x = x₁ and x = x₀.

Application of One Point Collocation Second Modification (mod 2) to a First-Order Reaction

Now we apply the rational function approximation to the diffusion reaction problem with first order reaction rate, where:

Let us choose C such that:

Thus,

We note that this rational function is exact for terms up to In addition, as u_r(x)→0 for x [0,1).

We note that the standard collocation has the same values for (u) at the collocation points as that of mod.1 and they differ at u(0). The mod. 1 and mod. 2 have the same values of (u) at the collocation points and at x = 0.

EFFECTIVENESS FACTOR CALCULATIONS

There are many ways to derive expressions for the effectiveness factor. Equation (10) is useful for the cases of low

If we multiply Equation (7) by and integrate both sides, we obtain:

For a slab, s=0;

This formula is useful for large value of . Going back to Equation (7) and multiplying both sides by the term

and then carry out an integration of both sides with respect to (r), we reach, after a lengthy derivation, the following expression for the effectiveness factor:

where

and

Solving the second order Equation (47) in , we obtain:

This simplifies for the case of a slab to:

and for a sphere to:

Note that I₁ equals zero for a slab and can be integrated analytically for a sphere. For a first-order reaction, I₁ and I₂ can be integrated analytically to give known expressions for the effectiveness factor.

NUMERICAL RESULTS

Example 1

The first example deals with a first-order reaction for which the application of the collocation method leads to a system of linear equations. For a first-order reaction occurring in a slab with = 100, we applied different methods discussed in this paper to plot the dimensionless concentration (u) against dimensionless distance (x). The results are shown in Figure 1. The one point standard collocation gives a negative concentration at x = 0, whereas the two point collocation profile oscillates. For the modified collocation methods, a very small negative value occurs at x = 0. The mod. 2 one point collocation profile is very close to the analytical solution.

The effectiveness factor is compared using different methods for different and is shown in Table 1. For one point standard collocation, the effectiveness factor is accurate using Equation (14) for small values of , whereas its accuracy improves as increases using Equation (46). The application of Equation (14) to the profile obtained by the application of mod. 2 to the one point collocation gives overall good accuracy, but will not be able to predict the asymptotic value of the effectiveness factor as increases.

The application of Equation (14) to the mod. 1 collocation method will give the same results as the one point standard collocation because the integral depends on the value of u at the collocation point, which is the same in the two cases.

The application of Equation (46) gives better accuracy and will converge to the asymptotic value for large . The results will be the same for mod. 1 and mod. 2 because the integral in Equation (46) depends on u(0), which is the same in both methods. The results for the two point standard collocation using Equation (14) are excellent up to = 20, but will start to deteriorate for higher values. Moreover, it can be noted from Table 1 that the results for the two-point standard collocation method using Equation (46) are better than those using Equation (14) for >10, but are slightly worse than those of mod. 2 one-point collocation for > 10.

In summary, the results show the superiority of mod. 2 one-point collocation for the calculation of the concentration profile and the effectiveness factor using Equation (46).

Example 2

We now consider a 4^th order reaction taking place in different shapes; slab, infinite cylinder and a

sphere. The results for the case of a slab are given in Table 2 for the effectiveness factor for different . Here Equations (46) and (51) are the same, and only terminal values of u are needed. One point collocation with mod. 1 or 2 and Equation (46) gives reasonable results, while exact results are obtained for two-point collocation with mod. 1 or 2.

For the case of a slab, it can be seen from Table 2 that for large Thiele modulus, the results for the standard two-point collocation method using Equation (46) take intermediate values between one and two point collocation with mod. 1 or 2 using Equation (46). The latter is the closest to the exact solution.

For the case of a cylinder (Table 3), the effectiveness factor calculations require an accurate concentration profile. The collocation with mod. 2 using Equation (51) gives the best results. The results for two-point collocation give, on the other hand, a reasonable accuracy. It can also be noted that substantial improvement is obtained with standard two-point collocation using Equation (14) over the one-point standard collocation. The same observations are noted for the case of a sphere in Table 4.

It can also be seen that, for both the cases of cylinder and sphere and for large Thiele modulus, the results for the two-point standard collocation method using Equation (51) give intermediate results between one and two point collocation with mod. 2 and using Equation (51). The latter is the closest to the exact solution.

The use of four-point collocation (not shown) gave exact results. Less accurate results are obtained when using Equation (45) for the evaluation of the effectiveness factor.

CONCLUSIONS

Two modifications are presented to improve the performance of the standard collocation method. In the first modification, we add an extra collocation point to the zeros of Jacobi polynomials at the center of the catalyst particle. The solution is in the form of a higher order polynomial. In the second modification this polynomial is transformed into a rational function with better accuracy of the solution. When this rational function is used with the proper formula for the effectiveness factor, excellent results are obtained. The results are of particular importance to the case of large Thiele modulus where the concentration profile is very steep. Such a case used to be treated by the dead-zone method (Soliman, 1989).

ACKNOWLEDGEMENT

The authors extend their appreciation to the Deanship of Scientific Research at King Saud University for funding the work through the research group project No RGP-VPP-188.

Submitted: July 5, 2013

Revised: January 27, 2014

Accepted: February 1, 2014

• Berrut, J. P., Baltensperger, R. and Mittelmann, H. D., Recent development in barycentric rational interpolation. In: Trends and Applications in Constructive Approximation. International Series of Numerical Mathematics, 15 27-51 (2005).
• Biscaia Junior, E. C., Mansur, M. B., Salum, A. and Castro, R. M. Z., A moving boundary problem and orthogonal collocation in solving a dynamic liquid surfactant membrane model including osmosis and breakage. Braz. J. Chem. Eng., 18, 163 (2001).
• Finlayson, B. A., The Method of Weighted Residuals and Variational Principles. Academic Press, New York (1972).
• Finlayson, B. A., Nonlinear Analysis in Chemical Engineering. McGraw-Hill, New York (1980).
• Soliman, M. A., A modified one point collocation method for diffusion reaction problems. Chemical Engineering Science 43, 1198 (1988).
• Soliman, M. A., Collocation with low order polynomials for fast reactions in catalyst particles. Chemical Engineering Science, 44, 1459 (1989).
• Torres, L. G., Martins, F. J. D. and Bogle, I. D. L., Comparison of a reduced order model for packed separation processes and a rigorous nonequilibrium stage model. Braz. J. Chem. Eng., 17, 915 (2000).
• Villadsen, J. V. and Stewart, W. E., Solution of boundary-value problems by orthogonal collocation. Chemical Engineering Science, 22, 1483 (1967).
• Villadsen, J. V. and Michelsen, M. L. Solution of Differential Equation Models by Polynomial Approximation. Prentice-Hall, New Jersey (1978).

*

To whom correspondence should be addressed

Publication Dates

• Publication in this collection
  14 Nov 2014
• Date of issue
  Dec 2014

History

• Received
  05 July 2013
• Accepted
  01 Feb 2014
• Reviewed
  27 Jan 2014