Simulation of Distributed Parameter Systems Using a Matlab-Based

Ind. Eng. Chem. Res. 2004, 43, 3469-3477

3469

Simulation of Distributed Parameter Systems Using a Matlab-Based Method of Lines Toolbox: Chemical Engineering Applications Alain Vande Wouwer,*,† P. Saucez,‡ and W. E. Schiesser§ Service d’Automatique, Faculte´ Polytechnique de Mons, 31 Boulevard Dolez, 7000 Mons, Belgium, Service de Mathe´ matique et Recherche Ope´ rationnelle, Faculte´ Polytechnique de Mons, Mons, Belgium, and Mathematics and Engineering, Lehigh University, Bethlehem, Pennsylvania 18015

Computational modeling is now generally accepted as an essential procedure for the dynamic analysis of chemical processes. Many of these processes are distributed parameter systems, i.e., systems in which state variables depend on several independent variables (such as time and space) and which are described by sets of nonlinear partial differential equations (PDEs). The method of lines (MOL) is probably the most widely used approach to the solution of evolutionary PDEs, and the objective of this paper is to report on the development of a Matlab-based MOL toolbox. The toolbox contains a set of linear spatial approximation techniques, e.g., finitedifference methods, implemented using the concept of differentiation matrices, as well as a set of nonlinear spatial approximations, e.g., flux limiters. In addition, several time integrators, including basic explicit methods and some advanced linearly implicit methods, are included. The underlying philosophy of these developments is to provide the user with a variety of easily understood methods and a collection of application examples that can be used as Matlab templates for the rapid prototyping of new dynamic simulation codes. In this paper, Burgers’ equation in one and two space dimensions, as well as a dynamic model of a three-zone tubular fixed-bed reactor used for studying benzene hydrogenation, and the poisoning kinetics of thiophene on a nickel catalyst are considered. 1. Introduction Partial differential equations (PDEs) play an important role in the mathematical modeling of engineering systems, particularly chemical engineering systems. PDE models, which are derived from mass, energy, and momentum balances, usually describe space and time variations of the state variables, that can occur for various reasons, e.g., nonideal mixing conditions in batch or batch-fed reactors, material transport, reaction, heat and mass transfer in packed columns, and tubular reactors. Even though space-time is the most common combination of independent variables, some other combinations might be used, e.g., individual particle size and time in dynamic population models used to describe crystallization or polymerization processes. One of the most popular methods for solving evolutionary PDEs is the method of lines (MOL), which proceeds in two basic steps (in the following, it is assumed that space and time are the two independent variables under consideration): (i) spatial derivatives are first approximated using finite-difference, -element, or -volume methods; (ii) the resulting system of semidiscrete (discrete in space and continuous in time) equations is integrated in time. The success of the MOL originates in the availability of efficient time integrators for solving a wide range of problems, including ordinary differential equations * To whom correspondence should be addressed. Tel.: +32-(0)65-374141. Fax: +32-(0)65-374136. E-mail: [email protected]. † Service d’Automatique, Faculte´ Polytechnique de Mons. ‡ Service de Mathe´matique et Recherche Ope´rationnelle, Faculte´ Polytechnique de Mons. § Lehigh University.

(ODEs) and mixed systems of algebraic equations (AEs) and ODEs (AEs and ODEs forming a system of differential-algebraic equations, DAEs). The MOL has attracted considerable attention, and consequently several general purpose libraries [e.g., the Fortran codes in DSS/21 and NAG2,3 and software packages dedicated specifically to chemical engineering applications (e.g., gPROMs4 and DIVA5)] have been developed that include several spatial discretization and time integration algorithms. Even though these libraries and software have proved very useful in many applications, there is still a need for basic, introductory, yet efficient, tools for the simulation of distributed parameter systems (DPSs), i.e., software tools that can be easily used by practicing engineers (whose objective is probably not process simulation as much as process analysis and improvement through timely simulation studies) and that would provide up-to-date numerical algorithms. The objective of this paper is to report on the development of a MOL toolbox within Matlab. This toolbox provides a series of spatial discretization methods, including finite-difference and -volume methods, spatial adaptive gridding, and a set of time integrators. Several examples, in specific application areas, can serve as templates for developing new programs. The idea here is that a simple code template is often more comprehensible and flexible than a software environment with specific user interfaces. In this way, various problems stated as coupled systems of AEs, ODEs, and PDEs in one or more spatial dimensions can easily be developed and tested. Specifically, the following applications and methods are discussed in this paper: (i) Burgers’ equation, which illustrates basic finite-difference and -volume methods;

10.1021/ie0302894 CCC: $27.50 © 2004 American Chemical Society Published on Web 09/18/2003

3470 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004

some results are also shown for the application of adaptive gridding techniques (i.e., spatial grid refinement); (ii) Burgers’ equation in two spatial dimensions, which demonstrates how easily the previous algorithms can be applied to problems in two or three dimensions (on simple geometries); (iii) benzene hydrogenation in a three-zone tubular reactor, experiencing catalyst activity decay due to poisoning by thiophene, which is a more complex problem involving several coupled mass and energy balance PDEs, a catalyst activity PDE, and sets of BCs and continuity conditions at the inert zone and catalyst interfaces. Matlab is now widely available in industry and academia and provides a very convenient basis for the development of a MOL toolbox (called Matmol), requiring minimum programming expertise. Most of the spatial approximation techniques, because they are linear operators, can be implemented using the concept of differentiation matrices; i.e., spatial derivatives of a dependent variable can easily be obtained through a matrix multiplication (which is the key operator in Matlab). However, not all of the spatial approximation techniques follow this format. Nonlinear approximation techniques, ensuring oscillation-free numerical solutions [total variation diminishing (TVD) schemes, such as slope and flux limiters], might be required to handle problems that develop sharp spatial gradients (e.g., temperature or concentration fronts). Both method classes are considered in the paper. Note that Matlab has recently been used as a basis for the development of high-quality libraries for the solution of ODEs and DAEs6 and for the solution of boundary-value problems (BVPs) and PDEs by spectral collocation methods.7,8 The paper is organized as follows. The next section introduces the basic principles of the MOL and addresses implementation issues within Matlab. Several finite-difference and -volume methods, as well as a static regridding algorithm (spatial grid refinement), that are representative of the current development status of the MOL toolbox are considered. Section 3 covers several application examples that illustrate several important concepts, including upwinding to handle strongly convective problems, spatial grid adaptation to capture sharp moving fronts, and BC formulation. Finally, section 4 is devoted to concluding remarks and perspectives for the development of the MOL toolbox. 2. Matmol: A Matlab-Based MOL Toolbox

0 ) b(t,z,x,xz) x(t)0,z) ) x0(z)

z ∈ D, t > 0 z ∈ S, t > 0 z ∈ D ∪ S, t ) 0

x˜ t ) f˜(t,x˜ )

(4)

0)b ˜ (t,x˜ )

(5)

x˜ (t)0) ) x˜ 0

(6)

where x˜ is a vector of discretized dependent variables (e.g,. the values of x on a spatial grid covering the spatial domain D ∪ S) and f˜ and b ˜ are vectors of nonlinear functions approximating the original equations (1) and (2), respectively. (2) Time Integration: the resulting system of semidiscrete (discrete in space and continuous in time) equations (4)-(6) can be integrated in time using an available ODE or DAE solver. In practice, many problems can be solved efficiently using relatively simple ODE integrators (this is one of the main reasons for the popularity of the MOL). In particular, the use of a more sophisticated DAE solver can be avoided by using the BC (5) to express some unknown variables that can be substituted back into eq 4. This analytical procedure is, however, limited to simple BCs, and more complex boundary value problems can require the use of a DAE solver, as we will see in one of the application examples considered in the next section. These two basic steps are now discussed with regards to their practical implementation within Matlab. 2.1. Spatial Discretization. Historically, finite differences have been the first choice approximations used to solve PDE problems. They are based on a set of spatial grid points zi, i ) 1, ..., n, distributed in the spatial domain (for the sake of simplicity, a onedimensional spatial domain and a uniform spatial grid are considered first) and Taylor series for several values x(zi+1), x(zi+2), ..., x(zi-1), x(zi-2), etc. Linear combinations of several of these Taylor series give specific formulas, e.g.,

|

Consider an initial BVP (IBVP) in the form

xt ) f(t,z,x,xz,xzz)

D, their associated boundary conditions (BCs) defined on the boundary surface S of D, and initial conditions (ICs) defined on the complete spatial domain. The MOL is a numerical procedure for solving IBVP (eqs 1-3) that proceeds in two separate steps: (1) Spatial Discretization: the spatial derivatives appearing in eqs 1 and 2 are replaced by numerical approximations, e.g., finite differences, elements, or volumes, in such a way that eqs 1-3 are transformed into a DAE system

(1)

x(zi+1) - x(zi-1) ∂x + O(∆z2) ) ∂z zi 2∆z with ∆z ) zi - zi-1 (7)

(2)

is a second-order (three-point) centered finite-difference approximation for the first derivative ∂x(zi)/∂z.

(3)

∂x ) ∂z zi -x(zi-3) + 6x(zi-2) - 18x(zi-1) + 10x(zi) + 3x(zi+1) + 12∆z O(∆z4) (8)

where x is the vector of dependent variables (e.g., temperature, concentration), z ) [z1, z2, z3]T is the vector of spatial independent variables (in general, three orthogonal spatial coordinates), and t is an initial independent variable (typically time). A subscript notation is used for the several partial derivatives, i.e., xt ) ∂x/∂t, xz ) ∂x/∂z, and xzz ) ∂2x/∂z2. Equations 1-3 represent a system of PDEs defined in a spatial domain

|

is a fourth-order (five-point) biased-upwind finite-difference approximation for the first derivative ∂x(zi)/∂z. This latter formula is particularly useful to approximate

Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3471

strongly convective terms, with a propagation from left to right (i.e., positive flow velocity v > 0).

|

∂ 2x ) ∂z2 zi -2x(zi-2) + 32x(zi-1) - 60x(zi) + 32x(zi+1) - 2x(zi+2) 4!∆z2

(ii) Direct differentiation: x˜ zz ) D2x˜ +

O(∆z4) (9)

is a fourth-order (five-point) centered finite-difference approximation for the second derivative ∂2x(zi)/∂z2. When using a uniform spatial grid, formulas such as eqs 7-9 can be used at all of the grid points, with the exception of a few points near the boundaries (the number of these boundary points depends on the approximation stencil of the formula, i.e., the number and configuration of the grid points involved in the formula), where alternative formulas have to be developed in order to avoid the introduction of fictitious grid points outside of the spatial domain. For instance, in the case of eq 7

| |

-3x(z1) + 4x(z2) - x(z3) ∂x + O(∆z2) ) ∂z z1 2∆z

(10)

3x(zN) - 4x(zN-1) + x(zN-2) ∂x + O(∆z2) (11) ) ∂z zN 2∆z which are no longer centered formulas. For the Matlab implementation of formulas such as eqs 7, 10, and 11, the concept of a differentiation matrix is useful, i.e.,

[

-3 4 -1 0 ... -1 0 1 0 ... ... 1 0 ... -1 0 1 x˜ z ) D1x˜ ) 2∆z ... 0 ... ... 0 -1 0 ... ... 0 1

][ ]

0 x˜ 1 0 x˜ 2 ... ... 0 x˜ i ... 0 1 x˜ N-1 -4 3 x˜ N

... ...

(i) Stagewise differentiation: x˜ z ) D1x˜ and x˜ zz ) D1x˜ z (13)

(12)

so that a spatial derivative can be computed by a simple matrix operation (which is the key operation in Matlab!). Similar differentiation matrices can be developed in the case of approximations (8) and (9), by using the appropriate boundary formulas. The matrix elements are the weighting coefficients involved in the several approximation formulas. These weighting coefficients can be computed analytically in a number of cases, i.e., various centered or noncentered computational stencils for first-, second-, and third-order derivatives; the weighting coefficients can also be generated automatically, in a numerical way, using an algorithm due to Fornberg,9,10 which gives the weighting coefficients of finite-difference formulas of arbitrary order of accuracy on arbitrarily spaced spatial grids (nonuniform grids). Both approaches have been included in Matmol. The first approach has a pedagogical value, whereas the second one focuses on generality and efficiency. Moreover, differentiation matrices allow higher-order spatial derivatives to be easily computed in either of two ways:

(14)

where D1 and D2 are matrices for computing first- and second-order derivatives, respectively. Other combinations are possible as well, e.g., x˜ zzzz ) D1{D1[D1(D1x˜ )]} or x˜ zzzz ) D2(D2x˜ ) or x˜ zzzz ) D4x˜ . Stagewise differentiation is a simple procedure, which can be very useful in some cases but which should be applied with care because an order of accuracy is lost at each stage. Many approximations, such as finite differences and finite elements, which are linear in nature, can be cast into the form of differentiation matrices. However, nonlinear approximation techniques may be required when studying PDE problems developing sharp spatial transitions (e.g., boundary layers, traveling waves) for which classical linear techniques either produce excessive numerical dispersion (in the case of first-order schemes) or spurious oscillations (in the case of higherorder schemes). To mitigate these problems, high-order oscillation-free schemes have been proposed for hyperbolic conservation laws of the form

∂f[x(t,z)] ∂x(t,z) )∂t ∂z

z ∈ D, t > 0

(15)

These schemes make use of a high-order method in smooth parts of the solution and modify it (in a nonlinear way) near sharp spatial transitions so as to increase the amount of numerical dispersion. These methods include flux- and slope-limiter schemes.11 In flux-limiter methods, the numerical approximation ˜f of the flux f in eq 15 can be viewed as the combination of a low-order approximation ˜flo and a correction

˜f(x˜ ) ) ˜flo(x˜ ) + φ(x˜ ) [f˜ho(x˜ ) - ˜flo(x˜ )]

(16)

where ˜fho is a high-order approximation of f and φ(x˜ ) is the limiter. The role of φx˜ ) is to sense the smoothness of x˜ and to accordingly switch between the low- and high-order approximations. Intuitively, if x˜ is smooth, then φ(x˜ ) should be near 1, whereas in nonsmooth regions, φ(x˜ ) should be near 0. In fact, however, a larger range of variation for φ(x˜ ) allows desirable properties to be achieved, as will be apparent in the subsequent discussion. There are various ways to measure the smoothness of the numerical solution, but a popular approach is based on the ratio of consecutive solution gradients, i.e., φ(ri) with

ri )

x˜ i - x˜ i-1 x˜ i+1 - x˜ i

(17)

If r is near 1, then the numerical solution x˜ is likely to be smooth, whereas smaller or larger values indicate the presence of spatial transitions. In work by Sweby,12 a variety of limiter functions are studied and conditions are derived that ensure that the flux-limiter method is second-order accurate and TVD. Loosely speaking, this latter property amounts to saying that the total variation of the numerical solution x˜ in both time and space is monotonically decreasing as time evolves (see work by LeVeque11 for a more rigorous description). The TVD property is ensured if


Figure 1. Several limiter functions ensuring second-order accuracy and TVD property.

0e

φ(r) e2 r

0 e φ(r) e 2

(18)

Second-order accuracy is obtained if φ(r) passes smoothly through the point φ(1) ) 1. In addition, Sweby12 recommends the limiter to be a convex combination of the Lax-Wendroff limiter φ(r) ) 1 and the Beam-Warming limiter φ(r) ) r, to avoid too much distortion of smooth solutions. Consequently, Figure 1 shows the region in which φ(r) must lie in order to give second-order TVD methods, as well as several wellknown candidates. The Minmod and the Superbee limiters delineate the lower and upper limits of the second-order TVD region, respectively. Matmol includes a variety of flux-limiter methods that can be used to solve PDEs possessing solutions with sharp spatial gradients. 2.2. Time Integration. Once the spatial derivatives have been replaced by well-chosen approximations, the resulting system of ODEs or DAEs in the form of eqs 4-6 can be solved using one of the high-quality solvers provided in the Matlab ODE suite.13 This suite of time integrators includes (a) explicit Runge-Kutta methods for nonstiff systems of ODEs in functions ode23 and ode45, (b) numerical differentiation formulas (NDFs), i.e., implicit multistep methods for solving systems of stiff ODEs or systems of DAEs of index 1, in function ode15s, and (c) a modified Rosenbrock method, i.e., a linearly implicit one-step method for solving stiff systems of ODEs, in ode23s. Stiff systems, i.e., systems with largely different time scales (e.g., chemical reaction systems with slow and fast kinetics), are best solved using an implicit method. These methods require the Jacobian matrix J ) ∂fi/∂xj, which can be computed by finite differences (using the function numjac) or by a user-supplied function. The Jacobian matrix can be a dense or sparse matrix, and in this latter case, the sparsity pattern can be specified in another user-supplied function. At this stage, the Jacobian matrix can be directly defined by the user as a sparse matrix of 0’s and 1’s, using, for instance, the function spdiags (which allows the diagonals of sparse matrices to be created). When the matrix structure is more complex and cannot be defined easily using spdiags, the sparsity pattern can be evaluated numerically through a single call to numjac with the specification of a threshold value that allows zero

elements (smaller than the threshold value) to be discriminated from nonzero elements (larger than the threshold value). The resulting matrix can then be converted into a matrix of 0’s and 1’s with the function spy. The solver ode15s allows not only systems of stiff ODEs but also systems of DAEs of index 1 [e.g., eqs 4-6, where AEs (5) result from spatial discretization of the BCs] to be solved. An original feature of ode15s is the solution of either an ODE system or a DAE system, and the difference is completely transparent to the user.14 In other words, ode15s recognizes automatically that the problem is a DAE system and computes consistent ICs [i.e., values for x(t)0) and xt(t)0) that satisfy the equations at t ) 0]. The available time integrators are easy to use in a MOL solution procedure, and because they all have a similar calling sequence, several solvers can easily be compared in terms of solution accuracy and computational costs. They also offer sophisticated features, such as event detection (zero-crossing of special functions called events). However, this sophistication results in relatively complex codes; e.g., ode23s has more than 700 lines. For pedagogical purposes, we provide a set of basic time integrators, e.g., Euler, modified Euler, leapfrog, or explicit Runge-Kutta methods. In addition, a few time integrators for systems of stiff ODEs, e.g., linearly implicit Runge-Kutta methods, are also included. 2.3. Adaptive Gridding. Even though most of the time integrators automatically adjust the time-step size (and possibly the order of the integration formula) in order to meet stability and accuracy requirements, the conventional MOL proceeds only in a semiautomatic way because the spatial grid points are held fixed for the entire course of the computation. For problems developing large spatial transitions, such as steep moving fronts, this conventional approach can be inefficient because a large number of uniformly distributed grid points is required to adequately capture the regions of high spatial activity. Unfortunately, most of the nodes are “wasted” in smooth parts of the solution, and it is therefore desirable to use procedures that adapt the spatial grid (that is, move or add/delete nodes) so as to concentrate them in the regions where they are needed. In this section, a simple static grid adaptation method is considered (for more details, see Vande Wouwer et al.15). The conventional two-step MOL procedure is extended in the following way: (1) approximation of the spatial derivatives on a fixed nonuniform grid, (2) integration of the resulting semidiscrete ODEs over a small period of time, (3) adaptation of the spatial grid, and (4) interpolation of the solution to produce ICs on the new grid; i.e., the first two classical steps now occur on a fixed nonuniform grid, time integration is then halted to adapt the number and location of the spatial grid points (step 3), and new ICs are computed by interpolation in order to restart time integration (step 4). The main advantage of this approach is that the PDE solution (steps 1 and 2) and the grid adaptation procedure (steps 3 and 4) are uncoupled. Hence, it is easily implemented within Matlab. The main disadvantages are (a) the time integration is halted periodically to adapt the spatial grid (resulting in a computational overhead due to the frequent solver restarts), (b) as the grid points are moving at discrete times only, they may


be ineffectively placed (resulting in large temporal gradients when a steep moving front crosses some of the grid points, and therefore the ODE solver is required to use extremely small step sizes to retain accuracy), and (c) an interpolation technique is required to transfer the data from the old to the new grid, thus introducing additional sources of numerical errors. Despite these drawbacks, experience has indicated that static regridding can be very useful to study PDE problems involving traveling waves and to resolve smallscale features with improved accuracy. The grid adaptation mechanism used in the static regridding function agereg is based on the equidistribution of a given monitor function m(x) subject to constraints on the grid regularity, i.e.,

mi-1(x) ∆zi-1 ) mi(x) ∆zi ) c, ∆zi ) zi+1 - zi (19) ∆zi 1 e eK K ∆zi-1

i ) 2, ..., N - 1

(20)

where mi(x) is a discrete approximant of the monitor function that can measure various solution features such as arc length or curvature, ∆zi is the local grid spacing, and c and K are two adjustable parameters affecting the grid density and regularity. Additionally, the accuracy of the spatial derivative approximations (e.g., using finite-difference techniques) and the stiffness of the semidiscrete system of differential equations are largely influenced by the regularity and spacing of the grid points, which indicates the importance of limiting grid distortion using spatial regularization procedures. In step 2, time integration of the semidiscrete system of stiff ODEs is preferably accomplished using a onestep implicit time integrator. The fact that the time integrator does not require the past history of the solution (as in the case of a multistep method) is a desirable feature in connection with the frequent solver restarts (a variable-order multistep method reduces the order to 1 at the time of reinitialization; a high-order one-step method retains the same order of accuracy during the reinitialization). In addition, an implicit method is required to handle large systems of stiff ODEs. 2.4. Matlab MOL Template. A collection of application examples is implemented following a simple, standardized, structure: (i) A main program, problem•main, where the PDE model parameters are set and defined as global variables (i.e., variables available in other program units), the differentiation matrices (in the case of linear spatial approximation techniques) are computed and are also defined as global variables, the parameters required for time integration are set and calls to an ODE/DAE solver are performed; output/plots are created at the end of the run. (ii) A function, problem•pdes, which is called by the ODE/DAE solver, and where the approximations of several spatial derivatives are computed (either by using the differentiation matrices initially evaluated in the main program or by repeatedly calling additional functions implementing nonlinear approximation techniques), and the PDE model is built. At this stage, it is particularly important to use a compact vector/matrix programming style in order to avoid unnecessary “for loops”, which are very time-consuming within Matlab.

Figure 2. Comparison of analytical (dotted lines) and MOL solutions (solid lines) at time t ) 0, 0.1, ..., 1 min using threepoint centered finite differences to approximate the convective term (n ) 201 and µ ) 0.001).

(iii) A library of general purpose functions, which implements various spatial approximation and time integration algorithms. These examples can easily be modified to implement new applications and test various methods. 3. Applications In this section, a few IBVPs demonstrate the use of the above-mentioned solution procedures and Matlab tools. 3.1. Burgers’ Equation in One Dimension. Consider Burgers’ equation in one dimension,1 i.e., a simplified momentum balance equation in the form

xt ) -(0.5x2)z + µxzz ) -xxz + µxzz

(21)

which is an interesting test problem for the following reasons: (i) It is nonlinear yet possesses several known analytical solutions that can be used to evaluate numerical solutions, e.g., x(t,z) ) 0.1e-0.05(z-0.5+4.95t)/µ + 0.5e-0.25(z-0.5+0.75t)/µ + e-0.5(z-0.375)/µ e-0.05(z-0.5+4.95t)/µ + e-0.25(z-0.5+0.75t)/µ + e-0.5(z-0.375)/µ

(22) (ii) The nonlinear convective term xxz produces interesting features, such as front sharpening, which can have a major effect on the performance of numerical methods. (iii) The kinematic viscosity µ can be adjusted so as to change the equation from predominantly convective (small µ for difficult problems) to predominantly dispersive (large µ for easy problems). Here, a relatively small value µ ) 0.001 is considered. In a first run, Burgers’ equation is solved using threepoint centered finite differences for the first derivative, five-point centered finite differences for the second derivative, and a spatial grid with n ) 201 uniformly distributed points. The semidiscrete system of ODEs is integrated in time using ode15s. The numerical solution is compared to the analytical solution (22) in Figure 2. Clearly, the use of centered finite differences for the


Figure 3. Comparison of analytical (dotted lines) and MOL solutions (solid lines) at time t ) 0, 0.1, ..., 1 min using five-point biased-upwind finite differences to approximate the convective term (n ) 201 and µ ) 0.001).

Figure 4. Comparison of analytical (dotted lines) and MOL solutions (solid lines) at time t ) 0, 0.1, ..., 1 min using Koren’s flux limiter to approximate the convective term (n ) 81 and µ ) 0.001).

convective term produces nonphysical oscillations. This undesirable effect can be alleviated through the use of upwind approximations, and Figure 3 shows the results obtained with five-point biased-upwind approximations for the first-order derivative. The numerical solution is in much closer agreement with the analytical solution. Alternatively, oscillation-free approximation schemes can be used, and Figure 4 shows the results obtained with Koren’s flux limiter and a spatial grid with n ) 81 points. Numerical results are very satisfactory but are obtained at the price of higher computational costs as a result of the repeated evaluation of the nonlinear flux limiter and a more difficult time-stepping process (a higher number of Jacobian evaluations is required, resulting in a higher number of function evaluations). The reduction of the order of accuracy (from 2 to 1) of the spatial approximation in regions of high solution gradients is apparent from Figure 4; i.e., the sharp corners of the front are somewhat smeared. The main advantage of these schemes is that they guarantee the absence of oscillations (which are still present in Figure 3, even though they are hardly visible). Finally, a numerical solution can be obtained using standard finite-difference approximations, e.g., five-

Figure 5. Superposition of the analytical MOL solution and MOL solution, which are indistinguishable, at time t ) 0, 0.05, ..., 1 min using an adaptive grid (bottom plot) and five-point biasedupwind finite differences to approximate the convective term (n j ) 94 and µ ) 0.001).

point biased-upwind finite differences, on adaptive grids. Figure 5 shows the almost perfect numerical solution plotted after each grid adaptation (here, the grid is adapted after fixed time intervals ∆t ) 0.05; note that it would be more efficient to adapt the grid after a fixed number of integration steps, but this is not so easily implemented with the Matlab solvers). The evolution of the adapted grid is also depicted in Figure 5. Grid points are initially concentrated in the two solution fronts. As time evolves, front sharpening occurs; i.e., the leftmost front, which corresponds to higher solution values, moves faster and eventually catches the rightmost front. The grid points then cluster in a single spatial region, following the movement of the steep front. An average number of grid points n j ) 94 is used. The computational load is about the same as that for the fixed-grid, finite-difference results with n ) 201 of Figure 3 (the reduction in the number of grid points by a factor of 2 is offset by periodic grid adaptation requiring the computation of new differentiation matrices and the reinitialization of the time integrator). 3.2. Burgers’ Equation in Two Dimensions. Consider Burgers’ equation in two dimensions1

xt ) -xxz1 + µxz1z1 - xxz2 + µxz2z2

(23)

which has the analytical solution given by

x(t,z1,z2) )

1 (z1+z2-t)/2µ

1+e

(24)

A relatively smooth case is considered here, corresponding to µ ) 0.1, and solved on a square spatial domain with n1 ) 21 and n2 ) 21, resulting in a total of 441 grid points. The spatial derivatives in eq 23 are very easily computed using the previously derived differentiation matrices (five-point biased-upwind finite differences for first derivatives and five-point centered finite differences for second derivatives), applied in each of the orthogonal directions z1 and z2. Figure 6 shows the numerical solution at time t ) 0.6, and Figure 7, the absolute value of the difference between the numerical and analytical solutions.


where cB and cT are the concentrations of benzene and thiophene, T is the temperature, and θA is the catalyst activity. In these expressions, v is the gas velocity, DB and DT are diffusion coefficients, λeff is the bed effective thermal conductivity, Fc and FG are the catalyst and gas densities, Fjcjp is the average volumetric heat capacity, is the bed void fraction, rB, rT, and rd are the rates of hydrogenation, thiophene chemisorption, and poisoning, -∆H is the heat of benzene hydrogenation, R is the wall heat-transfer coefficient, Tw is the wall temperature, and R and L are the reactor radius and length. The numerical values of all of these parameters can be found in work by Weng et al.16 Because the reactor is subdivided into three sections, eqs 25-27 have to be solved with rB ) 0 and rT ) 0 in the inert sections and the complete set of equations (25)-(28) has to be solved in the catalyst section. A total of 10 PDEs therefore have to be solved simultaneously. In z ) 0 and z ) L (reactor inlet and outlet), standard BCs apply

Figure 6. Numerical solution at time t ) 0.6 min.

cB(t,z)0) ) cB,in(t), cT(t,z)0) ) cT,in(t)

(29)

T(t,z)0) ) Tin(t)

(30)

cBz(t,z)L) ) 0, cTz(t,z)L) ) 0

(31)

Tz(t,z)L) ) 0

(32)

At the section interfaces, z ) L1 and z ) L1 + L2 [where L1 is the length of the entrance section and L2 is the length of the catalyst section; note that L1/L ) 0.28 and (L1 + L2)/L ) 0.47], special continuity BCs apply, e.g., Figure 7. Absolute value of the error in the numerical solution at time t ) 0.6 min.

3.3. Benzene Hydrogenation in a Tubular Reactor. Consider a tubular fixed-bed reactor consisting of three sections: an inert entrance section, an active catalyst section, and an inert outlet section.16 A dynamic model is developed to study benzene hydrogenation (an exothermic reaction) and the poisoning kinetics of thiophene on a nickel catalyst. This model includes material balance equations for benzene and thiophene, an energy balance, and an equation accounting for catalyst deactivation. Because hydrogen is in great excess in all of the experiments considered by Weng et al.,16 the hydrogen concentration is essentially constant.

∂cB ∂cB ∂2cB Fc ) -v + DB 2 + rB(θA,cB,T) ∂t ∂z ∂z

(25)

∂2cT Fc ∂cT ∂cT ) -v + DT 2 + rT(θA,cT,T) ∂t ∂z ∂z

(26)

FGcpG ∂T λeff ∂2T ∂T 2R (T - T) + ) -v + ∂t Fjcjp ∂z Fjcjp ∂2z RFjcjp w -∆H F r (θ ,c ,T) (27) Fjcjp c B A B ∂θA ) rd(θA,cT,T) ∂t

(28)

+ + cB(t,z)L1 ) ) cB(t,z)L1 ), cT(t,z)L1 ) ) cT(t,z)L1 ) (33) + T(t,z)L1 ) ) T(t,z)L1 )

(34)

+ + cBz(t,z)L1 ) ) cBz(t,z)L1 ), cTz(t,z)L1 ) ) cTz(t,z)L1 ) (35) + Tz(t,z)L1 ) ) Tz(t,z)L1 )

(36)

where L1 denotes the rightmost limit of the entrance section and L+ 1 denotes the leftmost limit of the catalyst section. Equations 33 and 34 express concentration and temperature continuity between sections, whereas eqs 35 and 36 express the continuity of mass and energy fluxes (here, diffusion and heat conduction is the same in both sections so that the diffusion and thermal conductivity coefficients do not appear). Similar expressions hold at z ) L1 + L2. If abrupt time variations are expected in the inlet (forcing) variables cB,in(t), cT,in(t), and Tin(t), then BCs (29) and (30) are best implemented as ODEs, e.g.,

dT 1 (t,z)0) ) [Tin(t) - T(t,z)0)] dt τ

(37)

where τ is a small time constant representing the dynamics of a temperature change at the reactor inlet. In fact, abrupt time variations are always the result of


Figure 8. Evolution of the temperature profile inside the reactor at t ) 0, 1, 2, ..., 40 min: (a) reactor preheating for 0 < t e 20 min; (b) benzene hydrogenation for 20 < t e 40 min.

Figure 9. Catalyst poisoning: evolution of the temperature profile inside the reactor at t ) 40, 50, ..., 190 min.

mathematical idealization, and τ models, in a very simple way, the actuator dynamics. The 15 other BCs and continuity conditions can be implemented as AEs, leading, after approximation of the spatial derivatives in the PDEs and BCs, to a large system of DAEs. Here, finite differences (five-point biased-upwind finite differences for first derivatives and five-point centered finite differences for second derivatives) are used with n1 ) 71, n2 ) 71, and n3 ) 21 grid points in the first, second, and third reactor sections, respectively. The resulting 560 DAEs can be solved efficiently using ode15s. Figure 8 shows the evolution of the temperature profile inside the reactor operated in two distinct modes: (a) reactor preheating with Tin(t) ) 160 °C (for 0 < t e 20 min) and (b) benzene hydrogenation with Tin(t) ) 160 °C and cB,in(t) corresponding to a mole fraction xB,in(t) ) 0.066 (for 20 < t e 40 min). During phase a, the temperature gradually decreases from the inlet to the outlet of the reactor, as a result of heat exchange with the reactor jacket. During phase b, an exothermic reaction takes place, leading to the formation of a “hot-spot” at about one-third of the reactor length. Figures 9-11 illustrate the effect of catalyst poisoning. From t ) 40 min, a small quantity of thiophene is fed to the reactor. Consequently, catalyst activity decays

Figure 10. Catalyst poisoning: evolution of the benzene concentration profile inside the reactor at t ) 40, 50, ..., 190 min.

Figure 11. Catalyst poisoning: evolution of the catalyst activity profile inside the central (active) reaction zone at t ) 40, 50, ..., 190 min.

(Figure 11), and traveling waves of concentration and temperature form (Figures 9 and 10). 4. Conclusions This paper reports on the development of a Matlabbased MOL toolbox. In its current status, the toolbox includes a library of finite-difference and finite-volume methods and several basic time integrators. Ongoing developments concern essentially nonoscillatory schemes, pseudospectral and spectral methods, finite-element methods, additional time integration, and adaptive gridding algorithms, as well as a collection of application examples in various fields. Two important features of the toolbox are (a) emphasis on method diversity, i.e., we expect that the user who has heard about a class of methods will find at least one representative algorithm in the toolbox that he/she can use and test on his/her PDE problem, and (b) emphasis on readability and understandability of the several Matlab functions, so that the toolbox can easily be used by practicing scientists and engineers. To this end, relatively simple, yet effective, methods are included in the toolbox, and codes are thoroughly commented. The toolbox is available upon request from the authors (Alain.VandeWouwer@ fpms.ac.be or [email protected]).


Literature Cited (1) Schiesser, W. E. The Numerical Method of Lines: Integration of Partial Differential Equations; Academic Press: San Diego, 1991. (2) Berzins, M. Developments in NAG Library Software for Parabolic Equations. In Scientific Software Systems; Mason, J. C., Ed.; Chapman and Hall: New York, 1989; p 59. (3) Pennington, S. V.; Berzins, M. New NAG Library Software for First-Order Systems of Time-Dependent PDEs. ACM Trans. Math. Software 1994, 20, 63. (4) Oh, M.; Pantelides, C. C. A Modelling and Simulation Language for Combined Lumped and Distributed Parameter Systems. Comput. Chem. Eng. 1996, 20, 633. (5) Ko¨hler, R.; Mohl, K.; Schramm, H.; Zeitz, M.; Kienle, A.; Mangold, M.; Stein, M.; Gilles, E. D. Method of Lines within the Simulation Environment DIVA for Chemical Processes. In Adaptive Method of Lines; Vande Wouwer, A., Saucez, P., Schiesser, W. E., Eds.; CRC Press: Boca Raton, FL, 1991; p 367. (6) Shampine, L. F.; Gladwell, I.; Thompson, S. Solving ODEs with Matlab; Cambridge University Press: Cambridge, U.K., 2003. (7) Weideman, J. A. C.; Reddy, S. C. A Matlab Differentiation Matrix Suite, ACM Trans. Math. Software 2000, 26, 465. (8) Trefethen, L. N. Spectral Methods in Matlab; SIAM: Philadelphia, PA, 2000.

(9) Fornberg, B. Generation of Finite Difference Formulas on Arbitrarily Spaced Grids. Math. Comput. 1988, 51, 699. (10) Fornberg, B. Calculation of Weights in Finite Difference Formulas. SIAM Rev. 1999, 40, 685. (11) LeVeque, R. J. Numerical Methods for Conservation Laws; Birkha¨user Verlag: Basel, Switzerland, 1992. (12) Sweby, P. High-Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws. SIAM J. Numer. Anal. 1984, 21, 995. (13) Shampine, L. F.; Reichelt, M. W. The Matlab ODE Suite. SIAM J. Sci. Comput. 1997, 18, 1. (14) Shampine, L. F.; Reichelt, M. W.; Kierzenka, J. A. Solving Index 1 DAEs in Matlab and Simulink, available at ftp:// ftp.mathworks.com/pub/doc/papers/dae. (15) Vande Wouwer, A., Saucez, P., Schiesser, W. E., Eds. Adaptive Method of Lines; CRC Press: Boca Raton, FL, 2001. (16) Weng, H. S.; Eigenberger, G.; Butt, J. B. Catalyst Poisoning and Fixed Bed Reactor Dynamics. Chem. Eng. Sci. 1975, 30, 1341.

Received for review April 7, 2003 Revised manuscript received July 14, 2003 Accepted July 22, 2003 IE0302894

Simulation of Distributed Parameter Systems Using a Matlab-Based

Recommend Documents