MWRtools: a library for weighted residual method calculations

Computers and Chemical Engineering 23 (1999) 1041 – 1061 www.elsevier.com/locate/compchemeng

MWRtools: a library for weighted residual method calculations Yi-hung Lin, Hsiao-Yung Chang, Raymond A. Adomaitis * Department of Chemical Engineering and Institute for Systems Research, Uni6ersity of Maryland, College Park, MD 20742, USA Received 29 May 1998; accepted 8 April 1999

Abstract We present a set of MATLAB functions developed for solving boundary-value problems using globally defined trial function expansions and weighted residual methods (MWR). The numerical techniques form a computational toolbox consisting of a common set of numerical tools for implementing the different MWR techniques used in the numerical solution and analysis of process systems described by partial differential equation models. This paper presents the results of our initial efforts to create computational modules that have a one-to-one correspondence with each step of implementing eigenfunction expansion, Galerkin projection, collocation, and other weighted residual methods. Examples illustrating the use of the library of subprograms are provided, including discussions of the development and application of more advanced MWR numerical techniques. © 1999 Elsevier Science Ltd. All rights reserved. Keywords: Boundary value problems; Weighted residual methods; Eigenfunction expansions; Galerkin projection method; Orthogonal collocation; Discretization methods

1. Introduction There has been a considerable amount of research devoted to developing numerical techniques for solving boundary-value problems (BVPs) in the chemical engineering literature (e.g. Villadsen & Stewart, 1967; Finlayson, 1972; Villadsen & Michelsen, 1978; Michelsen & Villadsen, 1981). The large size of this body of literature reflects the prevalence of BVPs in chemical process modeling problems involving transport and distributions defined on finite domains. For many of these problems, it is natural to pose solutions u(x, y, z, t) in the form of trial function expansions Nx,Ny ,Nz

u(x, y, z, t)=

%

ai, j,k (t)ci, j,k (x, y, z)

(1)

i, j,k = 1

where each trial function ci, j, k (x, y, z) is defined over the entire physical domain V and its boundaries #V. Partial differential equation solution procedures based on this form of trial function expansion fall under the broad classification of the spectral methods. In this paper, we focus on developing a set of computational * Corresponding author. Tel.: +1-301-4052969; fax: + 1-3013149920. E-mail address: [email protected] (R.A. Adomaitis)

techniques for implementing these discretization procedures. The residual functions to be minimized in the MWR procedure are found by substituting Eq. (1) into the partial differential modeling equation (u = G(u) (t

in V,

the boundary conditions F(u)= 0 that are not satisfied by the individual ci, j, k on #V, and the initial conditions u(x, y, z, 0)= u0(x, y, z) in V. There are a number of different methods for minimizing residual function norms and the most appropriate choice depends on the problem to be solved. The process of choosing the spatial mode amplitude coefficients ai, j, k to minimize the residual norm distinguishes the different weighted residual methods, but they all involve (exactly or approximately) the projection of the residual function R onto individual test functions fn using the weighted inner product R, fn w =

&

V

Rfnw(x, y, z)d6.

(2)

We note that this inner product definition separates the overall weight function into a general weight function

0098-1354/99/$ - see front matter © 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 9 9 ) 0 0 2 7 2 - 0

1042

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

w(x, y, z) and the weight associated with any assumed problem symmetries contained in the definition of the differential element d6. It is the method by which the projection operation is carried out or the definition of the test or weight functions that differentiates the MWR. Eigenfunction expansions are based on expressing the solution in terms of trial functions defined by the eigenfunctions of a Sturm–Liouville problem related to the BVP to be solved; the overall accuracy of this solution procedure normally depends on the accuracy to which the nonhomogeneous terms and initial conditions are approximated by the eigenfunction expansions. The Galerkin method is based on choosing a trial function expansion and projecting the residual onto each function, making the residual orthogonal to the sequence of trial functions. This method is applicable to nonlinear problems, especially when a pseudospectral formulation is used to simplify the implementation of the Galerkin projection. The least-squares projection is the optimal discretization technique for a given trial function expansion and is derived from minimizing the residual function norm. It is better suited to the solution of steady-state problems. Finally, the orthogonal collocation method is a pseudospectral technique based on a specific collocation point selection procedure that tends to reduce aliasing errors for nonlinear systems, giving a good and sometimes exact analog to the Galerkin method. The discrete-ordinate formulation of this method makes it convenient to use and provides much intuitive feel to the solution procedure.

1.1. Current MWR software There have been a large number of important contributions to the collection of numerical techniques available for solving BVPs by the different MWR. We present a partial survey of computational techniques based on globally-defined and spatially-localized trial function methods. Algorithms based on shooting techniques are not discussed. PDECOL is a package written by Madsen and Sincovec (1979) to solve systems of nonlinear partial differential equations in one space and one time dimension. The package implements finite element collocation methods based on piecewise polynomials for the spatial discretization. The boundary value problem solver COLSYS by Ascher, Christiansen and Russell (1981), or the enhanced version COLNEW, solves the mixedorder systems of BVPs defined by ordinary differential equations. In this spline collocation method, solutions are represented by B-spline functions evaluated at the Gaussian points defined in each subinterval. An adaptive spline point placement algorithm is used, based on minimizing a residual norm. The FORTRAN codes jcobi.f, intrp.f, dfopr.f and radau.f written

by Villadsen and Michelsen (1978) are a set of subroutines for solving BVPs with the orthogonal collocation method; the subroutines are used to compute the differentiation and quadrature arrays, and for interpolating solutions between the collocation points. The general chemical engineering software package Octave (Eaton, 1995) contains an orthogonal collocation subroutine to compute the collocation arrays based on the original algorithms of Villadsen and Michelsen (1978). The pseudo-spectral differentiation software package PseudoPack written by Don and Solomonoff (1997) computes differentiation arrays at the Gauss– Lobatto quadrature points using three different collocation methods. Fast transform and filtering algorithms are included in this package to enhance the speed and stability of computations. The MATLAB PDE toolbox MathWorks (1995) integrates a mesh generation algorithm, finite element technique, and a graphical user interface to solve BVPs in complex 1D or 2D geometries. The package is applicable to single partial differential equation models with constant coefficients and constant forcing terms. MacCluer (1994) includes a number of Mathematica examples of solving linear, one-dimensional BVPs, and Articolo (1998) presents a comprehensive overview of eigenfunction expansion solutions to one- and two-dimensional linear BVPs using Maple V. Analytical solution methods for linear and nonlinear partial differential equation problems are combined with MATLAB-based numerical techniques in Cooper (1998). A survey of BVP solvers that includes programs based on finite-difference, spline collocation, and projection methods with globally defined trial functions can be found in Childs, Scott, Daniel, Denman and Nelson (1979). Much emphasis has been placed on the numerical accuracy of computed quantities such as collocation array elements; as pointed out in Michelsen and Villadsen (1981), error analysis of computed solutions has been, for the most part, restricted to the comparison of numerical solutions to exact solutions. An example of one study that examined the convergence of solutions produced by BVP solvers is the paper of Denison and Hamrin (1983), which analyzed the numerical errors in BVP solvers DD04AD and COLSYS in the context of two test examples.

1.2. Organization of this paper This paper is organized in the following manner: an overview of our BVP solution framework is discussed in Section 2, together with the mathematical and computational elements that form the basis of our numerical discretization techniques; the MWRtools functions are presented in Section 3; eigenfunction expansion examples are discussed in Section 4; Galerkin projection examples are included in Section 5; and collocation

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1043

Fig. 1. The fine-scale discretization grid xˆ used to represent the trial functions. In this example, the trial functions are discretized equivalents of Bessel functions cn =J0(lnx) with J0(ln )= 0.

applications are discussed in Section 6. Higher dimensional spatial projection methods finally are discussed in Section 7.

function cn (x) defined over 05 x51 by the fine-scale trial function expansion M

cn (x)= % Cm,nlm (x)

(3)

m=1

2. The MWR solution framework The numerical techniques developed in this paper are distinguished from the studies cited in the previous section in that our goal was to develop a computational toolbox consisting of a common set of numerical tools for implementing the different MWR techniques inside the MATLAB computational environment. Our intention was to identify the numerical elements common to the different MWR and then to develop MATLAB functions that form a one-to-one correspondence between the subroutines and the elemental steps of an MWR solution procedure. Furthermore, we wanted to assess not only the accuracy of the elements of MWR implementation (such as the collocation arrays, inner product calculations, etc.) but also provide tools and techniques for assessing the accuracy of the solutions computed with these methods. The conceptual goal of making simple weighted inner product computations, computing eigenfunction sequences, etc., was to provide the next step in doing back-of-the-envelope type calculations for distributed parameter systems with methods for accurately assessing the discretization error.

2.1. Fine-scale discretization Because the trial functions (Eq. (1)) can be defined by polynomials, trigonometric functions, Bessel’s functions, or other functions, the key concept that made possible an interchangeable set of MWR computational tools was establishing a numerical technique for accurately representing these functions while retaining the maximum flexibility in choosing the form of the functions1. Consider, for example, representing the trial 1 Note that the same concept applies to the factors making up ci, j, k (x, y, z)=Xi (x)Yj (y)Zk (z) in higher-dimensional systems.

where the lm (x) are the building block polynomials of the Lagrange interpolation polynomial M

lm (x)=

(x−xˆi ) , (x ˆ i) i = 1, i " m ˆ m − x 5

m= 1, …, M.

(4)

The coefficient Cm, n represents the value of the function cn (x) at the discretization point xˆm because by definition (Eq. (4)), lm (xˆj )= 1 if m= j and lm (xˆj ) vanishes if m" j (see Fig. 1). If we now consider representing a state variable such as temperature u(x, t) using the truncated trial function expansion u N(x, t)=cb with c(x)= [c1(x), c2(x), …, cN (x)] b(t)= [b1(t), b2(t), …, bN (t)]T the discrete representation of the temperature approximation on the fine grid xˆ can be written as uM × 1 = CM × NbN × 1

with

Ci, j = cj (xˆi ).

Thus, the discrete transformation array C (the coefficients of interpolation function (Eq. (3))) is composed of columns corresponding to each trial function cn (x) evaluated at the fine-scale discretization points. For large M, it is numerically inappropriate to use directly the Lagrange interpolation polynomial to obtain values of u(x, t) between the fine-scale discretization points. Therefore, Neville’s algorithm (Stroud, 1974; Press, Flannery, Teukolsky & Vetterling 1988) is used to recursively construct the interpolated values.

2.1.1. Choosing xˆ and wˆ The method for choosing the number M and positions xˆm of fine-scale discretization points (Fig. 1) depends strongly on the problem and aspects of the MWR solution procedure. In general, the number M must be larger than the trial function truncation num-

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1044

ber N; for polynomial trial functions the minimum M is the maximum trial function degree q plus one, and for other trial functions the minimum M depends on the accuracy of the trial function approximation desired. Furthermore, accurate residual calculations for strongly nonlinear problems require more discretization points than are defined by the criteria above. For example, an accurate spectral decomposition of a residual generated by a quadratic function requires M = 2q + 1. We will demonstrate, however, that accurate inner product calculations do not necessarily require M = 2q + 1 (as a result of the product of two functions) by the choice of the quadrature weights. Many of the numerical difficulties encountered in solving BVPs are attributable to slow convergence rates of the residual in the neighborhood of the boundary condition locations. Thus, the choice of xˆm positions must balance interpolation performance against the accuracy of quadrature weights. Given all these considerations, we have chosen to fix the fine-scale discretization point locations as xˆ1 =0 and xˆM =1, and the interior points as the roots of the (M− 2)th-order a + 1) Jacobi polynomial J (1, (M − 2) (x), where a =0, 1, and 2, and corresponds to the slab, cylindrical, and spherical a + 1) geometries, respectively, with J (1, (x) = 1. These 0 polynomials are orthogonal with respect to weight function (1−x)x a + 1 in the unit interval. The roots are real and distinct, are all located inside the unit interval (Stroud, 1974), and are spaced in such a manner that the fine-grid discretization point density increases towards each end of the interval. This choice of fine-scale discretization points also guarantees that the quadrature weights wˆ, used to compute

&

1

a

T

f(x)x dx =w ˆ f.

0

where f. =[ f(xˆ1), …, f(xˆM )]T, result in exact integral evaluations when the degree of f is less than q = 2M− 3. Therefore, in the context of the inner product definition (Eq. (2)), the w ˆ m can be considered as the discrete approximations to the differential element d6.

2.1.2. Computing xˆ and wˆ Accurate and efficient algorithms for computing the values of xˆ as the roots of a Jacobi polynomial and the associated quadrature weights wˆ have been developed (Michelsen & Villadsen, 1972; Stroud, 1974). However, these algorithms can be modified to improve their numerical accuracy and efficiency, and to take advantage of vectorized MATLAB operations. Our algorithm starts by placing 3M points at the extrema of a (3M − 1)th degree Chebyshev polynomial, points that can be computed exactly and that have a distribution similar to the Jacobi polynomial roots we seek. Corresponding Jacobi polynomial values and their first derivatives are

evaluated with the general recursive formula that generates the Jacobi polynomials. The sign change intervals are identified and linear interpolation is used to determine initial guesses as the first step of root finding by the Newton–Raphson method. Roots that do not satisfy a specified accuracy tolerance are refined with the Newton–Raphson iterations until they do so. This defines the parallelized algorithm used to compute the root location vector xˆ. Since the xˆm are roots of a polynomial, the node spacing in the limit of the endpoints is inversely proportional to the square of the total number of points used. This suggests that these polynomial trial functions can resolve features with length scales of order M − 2 in the boundary layers while retaining good convergence properties in the bulk phase (see the discussions in Gottlieb and Orszag (1977), pp. 40). The Gauss–Lobatto quadrature weights w ˆ can be computed (Rice & Do, 1995) by w ˆ m=



CMKm dPM (xˆm ) dx



2

where PM (x)= M ˆ j ), Km = 1/(a + 1) for m= j = 1 (x−x 1, and Km = 1 for m"1. For a given M and geometry constant a, CM is computed as a function involving the product and quotient of Gamma functions2. However, as M increases, the machine accuracy of the operations used to compute CM soon saturate with round-off errors. To overcome this problem, we use the fact that M ˆ m = 1/(a + 1) must be true, thus, CM is comm=1 w puted from 1

CM = (a+ 1)%

M j=1



Kj dPM (xˆj ) dx



.

2

This algorithm eliminates the need to compute Gamma function values while improving the accuracy of the computed w ˆ m and extending the range of the maximum value of M.

2.1.3. Discrete differentiation and other operations Because the Lagrangian interpolation polynomials lm (x) are continuous and differentiable over the entire domain x[0, 1], explicit formulas for derivatives of functions represented by these trial function expansions up to order (M− 1) can be obtained. In the discrete-ordinate formulation of the collocation method, differentiation is a matrix operation, i.e. dT. /dx =A. T. and 92T. = B. T. where the elements A. j, i and B. j, i are defined as 2 The arguments to the Gamma function include M and a (see for example, Rice & Do, 1995).

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

A. j, i = and B. j, i =

dli (xˆj ) dx



dli (xˆj ) 1 d xˆ a xˆ aj dx j dx



for i, j=1, 2, …, M. The terms in each array are computed by directly differentiating the Lagrange interpolation equations, and the values of the nodal polynomial derivatives at the fine-scale discretization points xˆ are computed with the recurrence formulas derived by Michelsen and Villadsen (1972). The definition of the discrete Laplacian operator B. is then constructed from these differentiation arrays, using asymptotic values for the operator at the origin in the cases of cylindrical and spherical geometries to eliminate the numerical singularity caused by the 1/x a term.

2.2. Operations on the fine-scale grid The numerical discretization techniques previously described for constructing the fine-grid quadrature weights, differentiation arrays, and collocation point locations build on the computational techniques developed by Villadsen and Michelsen (1978). In our MATLAB implementation of the modified numerical techniques, we have found that these computations can be accurately carried out to approximately 500 discretization points. The effect of such finely discretized function representations is that numerical computations on this grid can be treated as (and under some circumstances are) exact, within the limits set by the number of discretization points M. Using this computationally accurate collocation method, we formulate our MWR solutions to BVPs in terms of a number N of trial functions that is smaller than the number of fine-scale discretization points M. This makes it possible to use nonpolynomial trial functions in the various MWR solution procedures, including collocation based on nonpolynomial trial functions. It also provides a means of estimating, or, for some problems computing exactly, discretization truncation errors and makes possible the development of exact discrete analogs to the Galerkin technique and other MWR methods (Adomaitis & Lin, 1998). We also note that some of the most computationally-intensive steps, such as Newton–Raphson iterations and the time integration of the discretized equations, depend primarily of the mode truncation number N and not the number M of fine-scale discretization points.

2.2.1. Elements of MWR By having an accurate method for representing arbitrary functions, and equally accurate methods for performing integration and differentiation operations, essentially all the MWR solution procedure operations

1045

reduce to matrix operations. For example, projecting the discretized function f. onto the N discretized orthonormal trial functions C is performed by aN × 1 = [CM × N]TDM × Mf. M × 1 where Di, j = 0 for i" j and Di, i = w ˆ i. These operations are performed by the weighted inner product function wip.m described in the following section. Other operations, such as the Gram–Schmidt orthonormalization of a sequence of functions and subprograms for computing collocation arrays from general sequences of orthogonal functions, all naturally build on the weighted inner product and other elemental operations. A time-consuming step in the implementation of an eigenfunction expansion solution procedure is the actual computation of the eigenvalues and definition of the eigenfunctions. In our solution framework, the sl.m function was developed (Adomaitis & Lin, 1999) to compute solutions (ln, cn ) to the subset of Sturm– Liouville problems defined by





1 d dc dc x ap(x) + q(x) +g(x)c= lc x 6(x) dx dx dx a

(5)

for x(0, 1). The solutions c(x) are subject to boundary conditions dc(0) dc(1) + b1c(1)= 0, a +bc(0)+a1 dx dx dc(1) dc(0) + dc(1)+ c0 + d0c(0)=0. dx dx

c

(6)

Therefore, Dirichlet, Neumann, mixed, periodic, and semiperiodic boundary conditions can be satisfied. The ln and cn are computed as the eigenvalues and eigenvectors of an array defined by the differentiation arrays corresponding to the differential operators and discretized functions vˆ, pˆ, qˆ, and gˆ. In the first step of this procedure, the boundary conditions are used to define a linear relationship between the values of c at the interval endpoints, and is used to reduce the discretized eigenvalue problem Eq. (5) to homogeneous form. After computing the M− 2 eigenvalues and eigenvectors, the endpoint values of the discretized eigenfunctions are computed to satisfy the boundary conditions, giving the matrix of discretized eigenfunctions CM × (M − 2). Error control is performed by evaluating the eigenvalue problem residuals generated when the eigenfunctions are interpolated to a finer discretization grid; eigenfunctions that produce residuals exceeding an infinity-norm bound are discarded (Adomaitis & Lin, 1999).

2.2.2. One-to-one correspondence to MWR To demonstrate a simple application of the numerical techniques developed, consider computing an eigenfunction expansion solution to the problem of deter-

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1046

mining the radial temperature distribution in a long, cylindrical rod with a uniform internal heat source and fixed temperature boundary condition T = 0 at r= 1. First, we pose the solution in terms of a trial function expansion N

T(r,t)= % an (t)cn (r) n=1

using the Laplacian operator eigenfunctions that vanish at the rod’s outer radial boundary and that are symmetric at r = 0. These eigenfunctions are the Bessel’s functions of the first kind of order zero, with eigenvalues determined from the zeros of this function. The residual formed by substituting the trial function expansion into the nonhomogeneous differential equation model

 

(T 1 ( (T = r +1 (t r (r (r

(7)

is minimized by projecting the residual onto each eigenfunction cn, giving explicit linear equations for the Bessel–Fourier coefficients of the expansion solution in time, subject to the initial condition T(r, 0) = 0. Each of these solution steps translates into a mathematical operation that involves solving a Sturm–Liouville problem, normalizing a set of trial functions, or projecting one function onto another with an inner product computation. As seen in Fig. 2, all of these operations can be performed a combination of the MWRtools functions and intrinsic MATLAB functions. An analogous comparison for the operations involved when implementing the Galerkin projection and orthogonal collocation methods can be developed for the same problem. As seen in Fig. 3, for both dis-

cretization methods we use the same set of polynomial trial functions cn (r) and test functions fn (r), where the latter are made orthonormal with respect to a weight selected to force the trial functions to zero at the outer boundary. Defining the residual functions and projecting this function onto the test functions appropriate to each method, the time derivative for the temperature profile is found in terms of the mode amplitude coefficients an or temperature measured at the collocation points Tcn. These sets of linear equations can be solved to find explicit solutions of temperature as a function of time. Five different MWRtools functions are used in the process of generating the trial functions and collocation arrays, performing the Gram–Schmidt orthogonalization, and computing the inner products needed for the projection operations in the above solution procedure.

3. The MWRtools functions The library of functions developed for implementing the weighted residual methods contains computational procedures for defining the problem domain and sets of trial functions, and general functions for spectrally discretizing problems, including collocation methods, that are described below. The notation describing the function input and output variables is Q for arrays, q for column vectors, q for scalar inputs, q for alphanumeric input, and {q} for parameters contained in MATLAB cell arrays.

Fig. 2. A comparison of the steps in an eigenfunction expansion solution procedure to Eq. (7) (left) compared to the corresponding computational steps of MWRtools (right). The solution is expressed in terms of a N = 5 mode eigenfunction expansion. M =30 points were used to define the fine grid used for the MWRtools solution.

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1047

Fig. 3. A comparison of the steps in Galerkin and collocation projection solution procedures to Eq. (7) (left) compared to the corresponding computational steps of MWRtools (right). The functions called to set up the physical domain and the method for reconstructing T(r, t) are the same as in the eigenfunction expansion example. Similarly, N = 5 and M= 30 for this example. 1 consists of a column vector of N 1’s.

[Q,w,A,B] = colmat(F,x,xˆ,wˆ,A. ,B. ) Computes the collocation differentiation arrays A and B, discrete transformation array Q, and quadrature weights w at the collocation points x based on the discretized set of trial functions F. Other inputs include the fine-grid discretization point locations xˆ, quadrature weights wˆ, and the discrete first-order differentiation and Laplacian operator arrays A. and B. , respectively. If only one fine-grid differentiation discretization array is specified in the input, the output is its corresponding collocation array. This function is normally called after colpts.m. [x] =colpts(F,xˆ,wˆ,xp,A. ) Computes collocation points x from a set of discretized trial functions F as the union of a set of pre-selected points xp and the roots of one of the trial functions. Roots are determined after the trial functions are orthonormalized with respect to quadrature weights wˆ if specified; the roots of each trial function are computed in a reverse sequence until the total number of collocation points equals the number of trial functions N. If xp is empty, N −1 roots are computed corresponding to an interior collection method. Since the roots are computed by the Newton – Raphson-based function rdf.m, increased accuracy can be obtained if the (optional) first-order differentiation array A. is provided.

[s] = fsf(n,N) Produces a vector of Fourier-space filter coefficients s used in the trial function expansions to reduce the oscillations associated with the Gibb’s phenomenon (Gottlieb & Shu, 1997). These oscillations result from projecting a discontinuous function (represented in discretized form) onto globally defined, smooth trial functions generated by either gdf.m or sl.m. The filter order is n, and N is the number of trial functions. [C] = gdf(xˆ,’f(xˆ,p)’,’xˆ ’,{’p’,p}) Generates discretized representations of a sequence of trial functions C according to the alphanumeric input formula f(xˆ,p), the fine-scale discretization grid xˆ and the parameter p. The cell array contents consist of the parameter name %p% and its numerical value(s) p. Additional parameters can be specified by concatenating the parameter name/value pairs. The functions can be polynomials, eigenfunctions generated by the explicit solution to a Sturm–Liouville problem, or an arbitrary sequence of functions chosen as part of a Galerkin discretization. This function is normally called after pd.m. [F] =gs(C,wˆ) A modified Gram–Schmidt orthogonalization procedure that provides better control of round-off errors than its classical counterpart. Starting with the first column (discretized function) in C, the orthogonalization algorithm proceeds by removing

1048

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

all components of the remaining functions not orthogonal to the current function in each iteration. Orthogonality is defined by the discretized weight function wˆ. Each trial function is also normalized during this process. Frequently, this function is used in conjunction with gdf.m. [xˆ,w ˆ ,A. ,B. ,Q. ] =pd(’geom’,M) Sets up the physical domain in terms of discrete points in the unit interval and defines the differentiation operators according to any imposed problem symmetries. This function is normally called first in the solution procedure. The inputs are the geometry (%slab%, %disk%, %sphe%, or %peri%) and the number of discretization points M. The output consists of the discretization grid xˆ, quadrature weights wˆ, the first-order differentiation array A. , and the discrete equivalent to the Laplacian operator B. . This function also computes the discretized set of M Jacobi polynomials Q. , a discrete transformation array used for filtering and spectral decomposition applications. The discretization point locations and quadrature weights are based on the modified Gaussian quadrature method (Gauss – Lobatto method). The two end points (0 and 1) are preassigned and interior positions are the roots of the (M−2)th degree Jacobi + 1) polynomial J (1,a M − 2 (x) where a is the geometry factor. The quadrature weights are computed after the discretization grid point locations are determined. The differentiation arrays are based on the Lagrange interpolation polynomial with the given M points. In this study, the recursive formula of the polynomial is used to construct the differentiation arrays. [y,dy] =polyint(xi,yi,x) Lagrange polynomial interpolation based on Neville’s algorithm. This recursive formula is obtained by rearranging the order in which the implementation of Lagrange interpolation calculations are performed. The result is spectrally accurate and provides an error estimate dy. The function takes the input interpolants yi with respect to xi and the positions x to be interpolated, producing the interpolated result y. [r] = rdf(x,f,df,e) Computes the roots r of a discretized function f. This function can be used to compute eigenvalues described by nonlinear relationships, and is also used for computing collocation point locations in colmat.m. The inputs consist of any discretization grid x and the corresponding function values f (and the optional first order derivative values df). The intervals bracketing the roots are found first and a linear interpolation in each interval is performed to find the roots. If the first-order derivatives are given, the roots are refined using the Newton – Raphson iterations (and polyint.m) until given tolerance e is satisfied.

[l,C,F,wef,wad] = sl(’geom’,A. ,xˆ,a,b,c,d,wˆ,vˆ,pˆ,qˆ,gˆ,errtol,a1,b1,c0,d0) A Sturm–Liouville problem solver, for problems of the type described by Eqs. (5) and (6), that computes the vector of eigenvalues l, discretized, orthonormal eigenfunctions C, adjoint eigenfunctions F, discretized weight function wef used to define the orthogonality of the eigenfunctions, and discretized weight function wad corresponding to the adjoint eigenfunctions. Input includes the problem geometry factor %geom% (geom= %slab% or %peri% corresponds to a= 0; geom = %disk% to a=1; and geom = %sphe% to a= 2), and the differentiation, discretization point, and quadrature arrays, A. , xˆ, and w ˆ , respectively. The remaining parameters correspond directly to the terms in Eqs. (5) and (6) or should be otherwise self-explanatory (Adomaitis & Lin, 1999). This function is normally called after pd.m. The maximum value errtol of the infinity-norm bound is used to discard inaccurate eigenfunctions; a negative value disables error control. ˆ) [Ip]= wip(F. ,G. ,w The weighted inner product of two sets of discretized functions F. and G. . This is one of the most heavily used subprograms since inner products are used in nearly every MWR step. The array of inner products Ip is computed as the vector dot product of the quadrature weight array w ˆ with the term-by-term (Hadamard-Schur) product of all combinations of the discretized functions represented in F. and G. . The number of rows in Ip corresponds to the number of trial functions in G. ; the columns correspond to the functions in F. .

4. Examples of Eigenfunction expansions Eigenfunction expansion solution procedures can be used for linear, nonhomogeneous problems where the nonhomogeneous terms appear in the modeling equation or the boundary conditions. Choosing trial functions as the eigenfunctions of the linear operator decouples the modes, thus, solutions can be computed term-by-term until a specified residual error tolerance is satisfied. Furthermore, the residual decreases monotonically with increasing truncation number.

4.1. Catalyst pellet poisoning Consider the spherical catalyst pellet with nonuniform catalytic activity (Carberry, 1976). By assuming the rate at which the poisoned section moves is slow relative to the rates of diffusion and reaction, we have the quasi steady-state problem

 

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1 ( 2(c r − k(r,t)f 2c =0 r 2 (r (r

(8)

subject to boundary condition c =1 at the pellet outer boundary (r =1) and #c/#r =0 at the pellet center (r=0). The slowly-moving catalyst activity profile is modeled using the catalyst activity function 0BkB 1, empirically modeled as a function of time with

1049

an = k(r,t)f 2,cn /ln which is computed as a=wip(phi^2*k,efun,w)./eval; conc=1+ efun*a; R=B*conc-k*phi^2.*conc;

errtol=20;

Reactant concentration profile snapshots for t=1, 10, 18, and 30 days are given in Fig. 5, showing the overall loss of reactivity of the catalyst pellet as a function of time. The solution residual function R N(r, t) is computed by substituting the truncated trial function expansion solution into the modeling Eq. (8). The residual L 2 norm is then computed as sqrt(wip(R, R, w)) and plotted as a function of N in Fig. 5. We observe that the norm decreases monotonically with increasing trial function expansion truncation number N. This demonstrates the concentration profile snapshot is a convergent solution, and can be used to determine the number of trial functions necessary to satisfy a desired error tolerance. The residual analysis also indicates that the error tolerance used to compute the discretized eigenfunctions does not adversely affect the solution convergence.

[eval efun] =sl(’sphe’,A,r,1,0, . . .

4.2. A distributed parameter system control problem

0,1,w,1,1,0,-phi^2*k,errtol);

A standard model for demonstrating the modal control of distributed parameter systems is the metal slab heating problem of Ray (1981), Ogunnaike and Ray (1994), and Gay and Ray (1995)

k(r, t)= 1/2 − tan − 1(200(r +t/20 −1))/p (see Fig. 4). The Thiele modulus for this problem is f= 20. We take solutions of the form c= 1 +N i = 1 aici where the ci are computed from the solutions to the eigenvalue problem lc =92c −k(r, t)f 2c subject to #c/#r =0 at r = 0 and c =0 at r =1. The eigenvalues and eigenfunctions are numerically computed by [r,w,A,B] =pd(’sphe’,100); t = 10; k = 0.5-atan(200*(r + t/20-1))/pi;

for the case of t= 10 days. We note that a large error tolerance is used on account of the numerical difficulties introduced by the sharp poisoning profile. Substituting the trial function expansion into the modeling equation to form the residual and then projecting the residual onto each normalized eigenfunction gives the solution as N

% lnancn =k(r,t)f 2 n=1

(y ( 2y (y = 2 − V − by(x,t)+ u(x,t) (t (x (x which describes heat transfer by conduction and a convective contribution due to the movement of the metal slab (see Fig. 6). With the representative parameter values V=10 and b= 20 taken from Gay and Ray

Fig. 4. The catalyst pellet (left) and the reactant rate activity coefficient k(r, t) as a function of pellet position and time, shown for t= 1, 10, and 18 days (right).

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1050

Fig. 5. Catalyst pellet reactant concentration profile snapshots (left) for t=1, 10, 18, and 30 days. Residual L 2 norm as a function of truncation number N for t = 10 days (right).

(1995), we assume solutions take the form y N(x, t)= fa and that the spatially-varying manipulated variable used for slab temperature control can be represented as u N(x, t)=fb. The trial functions fn are computed as the eigenfunctions of the non-self-adjoint eigenvalue problem lf =f%% − Vf%− bf subject to f% =0 at x= 0, 1. Defining the residual N

N

N

n=1

n=1

n=1

% a; nfn = % lnanfn + % bnfn and projecting this function onto the normalized adjoint eigenfunctions cn gives a; n = lnan +bn,

n= 1, …, N.

(9)

In this system, the control problem is defined as determining the values of u1(t), u2(t), and u3(t) of the uniformly-spaced heater inputs so that the slab temperature profile is maintained as nearly to yspt =(1− cos(px))/2 as possible. The modal control technique from Ray (1981) is used to determine u(x, t). The normalized eigenfunctions and eigenvalues are computed by [x w A] =pd(’slab’,50);

ing wip(e, c, w). The modes are fed into a PI controller to determine the control input function u(x, t) which is translated into corresponding control element commands ui, i =1,2,3 by projection. The synthesized control profile u(x, t) is then fed back to the process using wip(u, c, w) to determine the true bn used in the process dynamics simulations. Eq. (9) is numerically integrated to obtain the process response profile y(x, t)=fa. The closed-loop response of the system is shown in Fig. 8 for the case of N= 9. The PI controller gain of 0.5 and time constant 0.02 used in Gay and Ray (1995) were adopted for this study. The L 2 norm of e can be computed using the wip.m function to show monotonically decreasing behavior of e with time. This indicates the heating profile asymptotically approaches the optimal profile that is feasible with the discrete heating element inputs. Finally, we note that the singular function discretization methods discussed in Gay and Ray (1995) can be implemented using the MWRtools functions.

5. Galerkin projection examples

V = -10; beta=-20; [eval efun efun – ad] = sl(’slab’, . . . A,x,1,0,1,0,w,1,1,V,beta); In this case, 22 eigenfunctions were generated by sl.m that satisfy the default error tolerance. The first three eigenfunctions and adjoint eigenfunctions are plotted in Fig. 7. Simulation of this system under feedback control makes use a number of the MWRtools functions. The control loop consists of a comparator for computing the error e(x, t)=yspt(x) − y(x, t), where we assume that there are a sufficient number of temperature measurements to measure y(x, t) accurately. The modal analyzer projects e onto the adjoint eigenfunctions us-

When using the Galerkin discretization procedure, one represents the solution in terms of an orthogonal function expansion where each trial function satisfies the boundary conditions. The residual is projected onto these trial (test) functions in a process similar to the

Fig. 6. The three-zone moving metal slab temperature control problem.

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1051

Fig. 7. Eigenfunctions (left) and adjoint eigenfunctions (right) computed for the slab heating problem (cf., Fig. 11 of Gay and Ray (1995))

eigenfunction expansion method. However, this discretization procedure does not necessarily produce a modal decomposition that generates a decoupled set of ordinary differential equations in time; in fact, the mode amplitudes may be coupled in a nonlinear manner for nonlinear problems. This discretization method produces a convergent solution when the residual function can be represented by a convergent expansion in the trial functions cn. The tau method (Gottlieb & Orszag, 1977) is used when the trial functions do not automatically satisfy the boundary conditions. Numerical implementation of the semidiscrete Galerkin projection discretization procedure for nonlinear systems is typically carried out inside a function called by a MATLAB ODE solver. This function should be set up to reconstruct the state variables on the fine computational grid at each point in time and then use this representation of the states to evaluate the time derivative functions (in spatial position). The derivative functions are then projected onto the test functions using wip.m to find the corresponding mode amplitude coefficient time derivatives. Examples demonstrating the implementation of the Galerkin projection for a linear and nonlinear system follow.

5.1. Entry length problem Consider the entry-length problem described in Finlayson (1972) where a soluble gas species is absorbed by a falling liquid film (see Fig. 9). The absorbed species dimensionless concentration profile y(x, z) is described by the modeling equation 3 (y ( 2y − (1−x 2) = 2 2 (z (x with boundary conditions

(10)

Note the jump discontinuity in the boundary conditions at the corner point x= 0, z= 0. This discontinuity signals potential problems in solution convergence, therefore, the convergence of computed solutions in the neighborhood of this point must be checked. We take solutions in the form of the truncated trial function expansion N

y(x,z)= % an (z)cn (x)

(11)

n=1

where the individual trial functions cn (x) satisfy the homogeneous boundary conditions at x=0, 1. Using the Gram–Schmidt orthogonalization, we make the cn orthonormal with respect to weight (1− x 2). These operations are performed with the MWRtools functions nd – x=90; [x w A B] =pd(’slab’,nd – x); [eval efun] =sl(’slab’,A,x,0,1,1,0,w); psi=gs(efun,w.*(1-x.^2)); Substituting the trial function expansion (Eq. (11)) into the modeling Eq. (10) and projecting the residual onto each trial function using the inner product without the weight function (1− x 2) gives, a; = Fa with Fi, j = −

#

2 ( 2cj ,c 3 (x 2 i

$

(the dot denotes differentiation with respect to z) subject to an (0)= Ž1,cn )w computed with the inner product weight function. F can be factored into ULUT since F is symmetric, and the solution can be written directly as y(x,z)=cUe azUTa0

with a0n = an (0).

y = 0 at x = 0

The concentration profile at a representative value of z (e.g. z= − 0.1) can be computed by the following steps

(y/(x = 0 at x =1

N=40; psi=psi(:,1:N);

y= 1 at z= 0.

F=(-2/3)*wip(B*psi,psi,w);

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1052

a0 = wip(ones(size(x)),psi,w.*(1-x.^2)); z = -0.1; [U,Lam] = eig(F); y = psi*U*diag(exp(diag(Lam)*z))*U’*a0; A solution corresponding to a 40-mode trial function expansion is plotted in Fig. 9.

5.1.1. Solution con6ergence at z= 0 Taking a closer look at the contour plot in Fig. 9, we can observe oscillations in the computed solution near z = 0. This is not an indication of numerical problems, but is in fact an indicator of the Gibbs phenomenon (Gottlieb & Orszag, 1977), which is due to the discontinuity of the two boundary conditions at (x, z) = (0, 0). Using the polyint.m and wip.m functions of MWRtools, we present the results of convergence tests in Figs. 10 and 11. The left plot of Fig. 10 shows successive, norm-convergent approximations to the inlet condition by the projections onto sets of trial functions that vanish at x = 0. These curves are interplolated from the fine scale grid point values using polyint.m. The right plot of Fig. 10 demonstrates the pointwise convergence of the solution at z =0. Three points not on the fine grid are selected and polyint.m is used to compute the values at the selected points for a given N. As N increases, each series converges with faster convergence rates for points farther from x = 0. Finally, Fig. 11 (left) shows the monotonic convergence of the inlet condition residual in a normed sense computed with wip.m. Pointwise convergence performance can be improved using spectral filtering methods (Gottlieb & Shu, 1997). The coefficients sn used in the partial sum N

y(x,z)= % snan (z)cn (x) n=1

are computed using fsf.m: s = fsf(2,N);

yf=psi*U*diag(exp(diag(Lam)*z)) … *U’*(a0.*s); The visibly smoother results are presented in Fig. 11 (right); these curves correspond to the concentration profile at z= 0.

5.2. Compressor dynamics Axial compression systems are used in chemical process industry gas compression applications requiring large volume throughputs with relatively low pressure rise. Fluid catalytic cracking units are among the largest axial compressor users: air service for catalyst fluidization requires capacities of 50 000–300 000 cfm with pressure ratios of 2–4 over the atmospheric feed pressure. Other chemical and bulk material processing axial compressor applications include providing feed air for butadiene production by oxidative dedydrogenation of n-butene, nitric acid production by ammonia oxidation, air-separation plants, blast furnaces, and co-generation combustion applications. A schematic diagram illustrating the cross-section of a multistage axial compressor is shown in Fig. 12. Air enters at atmospheric pressure though the inlet diffuser and is compressed by a sequence of stages consisting of alternating rows of stationary and rotating blades. Flow control is provided by a combination of rotor speed adjustments, angle of the first rows of stator and inlet guide vanes, and throttle valve opening (located downstream of the exit diffuser). Overall pressure rise is determined by downstream process requirements and throttle valve pressure drop. The model of axial velocity component gas flow dynamics discussed in Adomaitis and Abed (1996) is used to describe the growth of small flow field perturbations to fully-developed stall cells after the stall bifurcation point (when the system is operated such that the intersection of the compressor characteristics falls to

Fig. 8. DPS control of a moving metal slab heating problem. (Left) closed-loop response to setpoint profile yspt(x)= (1 −cos(px))/2 (thick curve). The dashed curves represent slab temperature profiles from time =0 – 5. The control actions (right) are depicted as dashed curves and represent computed controller commands; the solid curves result from projecting the control action onto the three zone heater arrangement.

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1053

Fig. 9. Absorption of a gas into a falling film (left, after Fig. 3.1 of Finlayson (1972)). Dimensionless concentration of absorbed species for the entry length problem (right). Note the appearance of the Gibbs phenomenon in the upper right solution contour plot.

the left of the stall point). The model describes the dynamics of the annulus-averaged (mean) axial flowrate V(t) and the flow perturbation term 6(h, u, t) measured at the compressor face h= 0: Dp = f(V+60)−lc



dV ( − md dt (t

n

&

0

−

6dh

1 (6 (6 ( 26 − a 2 0 + 0 − m 20 . 2 (t (u (u In this model, 60 is the axial velocity perturbation evaluated at h = 0 (the inlet face of the compressor), Dp is the discharge-to-atmosphere pressure rise, h and u are the axial and angular coordinates, respectively, m = 0.001 is taken to be a gas eddy viscosity, a = 1/3.5 is the internal compressor lag, lc =8.0 is the overall compressor length, and md =1.75 is the exit diffuser length factor. Every compressor has its individual characteristic pressure rise f(Vloc) (Fig. 12). This function is obtained from experiments performed in the stable operating range and is estimated in the nonuniform-flow range. Following Moore and Greitzer (1986), we use a cubic equation in axial velocity

 

f(Vloc)=f0 + h 1+

 

3 Vloc 1 Vloc −1 − −1 2 w 2 w

n 3

where Vloc = V+60, the total local axial flow at the compressor face. The characteristic parameters are h= 0.18 as the pressure rise scaling factor, w=0.25 is the mean velocity scaling factor, and f0 = 0.3 is the shut-off head. A material balance on the gas in the diffuser and plenum chamber gives V= g Dp when the combined volume of the discharge diffuser and plenum is small relative to the compressor volume (this prevents surging flow (Moore & Greitzer, 1986)); this function represents the throttle characteristic and describes the overall resistance to flow generated by the compressor and throttle valve (see Fig. 12). The parameter g is proportional to the throttle opening. Following Moore and Greitzer (1986), if the flow disturbance is assumed to be irrotational it must satisfy Laplace’s equation, therefore, the overall trial function expansion for the flow field can be written as N

V(t)+ 6(h,u,t)= % an (t)e − ln hcn (u)

(12)

n=1

where the cn (u) are computed as nontrivial solutions to the eigenvalue problem d 2c du 2 subject to periodic boundary conditions. These trial functions and eigenfunctions can be conveniently comlc=

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1054

puted by using the combination of pd.m and sl.m functions: N = 2*10+1; nd – th = 6*N + 1; [th,w,dth,ddth] =pd(’peri’,nd – th); [lambda,psi] =sl(’peri’,dth,th ,… 0,1,1,0,w,1,1,0,0[ ],0,-1,-1,0); Defining V+ 60(u, t) =a1 +fa, the residual can be written as lc

da1 1 + fa; + mdfHa; dt a

=f(a1 + fa)−



n

d2 a 21 1 d −m 2 fa − 2 g du 2a du

(13)

with

adot(2:end) =H¯adot(2:end); In this function, the right side of Eq. (13) is evaluated on the spatial discretization grid using the velocity profile V+ 60(u) at the current point in time; this discretized function is projected onto the orthonormal functions cn using wip.m. The mode amplitude coefficients are then adjusted according to the left side of Eq. (13). Representative simulation results are presented in Fig. 13. This simulation begins with a small spatial perturbation that grows to a fully developed rotating stall cell. In Fig. 13, the growth is depicted in terms of the mode amplitude dynamics associated with the mean and first pair of spatially varying modes; snapshots of the fully developed stall cell are also included, illustrating the traveling wave nature of this flow instability.

H = diag(md./sqrt(-lambda(2:N)) + 1/alpha); subject to the constants (md and a) defined in the model description. Having computed the trial functions, quadrature weights, and model parameter values, the semi-discretized model is integrated forward in time using the ODE solvers supplied by MATLAB. The entire function containing the Galerkin projection is V = a(1); v0 = psi*a-V; rhs =f0+h*(1+ (3/2)*((V + v0)/omega-1))… -((V +v0)/omega-1).^3)/2 … -V^2/gamma^2 … -(dth-mu*ddth)*v0/(2*alpha); adot=wip(rhs,psi,w); adot(1) =adot(1)/lc;

6. Collocation examples The orthogonal collocation method (Villadsen & Stewart, 1967) was developed originally as a stable, predictable, and simple to implement pseudospectral technique. Because of its reliability, it has become a standard method for solving boundary-value problems by polynomial trial function expansions (Finlayson, 1972, 1980; Villadsen & Michelsen, 1978; Michelsen & Villadsen, 1981; Rice & Do, 1995). The interior formulation of this method (Villadsen & Stewart, 1967) is based on choosing a set of trial functions {cn }N n = 1 from an orthogonal polynomial sequence, with the discretization points computed as the roots of the polynomial cN + 1 next in the sequence. In most applications, orthogonal collocation approximates the Galerkin procedure, because the residual is forced to have as its

Fig. 10. Convergence analysis of the solution to the entry length problem at z =0. Solution curves y(x, z= 0) for N = 40, 50, and 60 are plotted on the left showing the Gibbs phenomenon near x = 0. The marker points indicate the solution at the computational grid points; interpolated values are obtained using polyint.m. (Right) pointwise convergence analysis for x = 0.05, 0.1, and 0.5.

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1055

Fig. 11. Inlet condition (z= 0) residual L 2 norm computed as a function of truncation number N (left). A comparison of original and spectrally filtered solution profiles at z =0 obtained using a 2nd-order Fourier-space filter (right).

primary component the polynomial used to determine the collocation points. However, for some linear problems where the residual can be expressed exactly in terms of the chosen set of trial functions, the two methods give identical results. The mixed collocation formulation is commonly used for problems with nonhomogeneous and nonlinear boundary conditions, and for BVPs defined by a set of differential equations with different boundary conditions. As in the case of the interior collocation formulation, the residual at N interior points is forced to vanish making the residual defined in V approximately orthogonal to the first N trial functions; additional collocation points are placed at the boundaries to satisfy the p boundary condition residuals, resulting in an overall trial function truncation number of N+ p. The result is that the mixed collocation method is a discrete approximation (and in some cases exact analog) to the tau method (Gottlieb & Orszag, 1977).

6.1. Orthogonal collocation arrays In the original orthogonal collocation reference (Villadsen & Stewart, 1967), the authors hand-calculate three example tables for collocation point locations, quadrature weights, and first-order differentiation and Laplacian operator arrays (see Table 2 of Villadsen and Stewart (1967)). The table entry for a given N (the number of interior collocation points) can be generated using four of the MWRtools functions: [x w dx ddx] =pd(’slab’,20);

array, and the collocation point locations for this mixed collocation example. The trial function orthogonality used to compute the collocation points includes the weight function (1− x 2). The accuracy of the discretization arrays computed with the MWRtools functions was found to be in the range of machine round-off error, matching the values given in Villadsen and Stewart (1967) to reported accuracy. We note that in this example, the trial functions used to generate Q are not multiplied by the function (1− x 2); in the following example, we consider the case where the trial functions contain this factor. With minor modifications, similar commands can be used to generate the Tables 5.2–5.5 in Finlayson (1972) for the orthogonal collocation method; the same set of MWRtools can be used for the least-squares-type collocation described in Theodoropoulou, Adomaitis and Zafiriou (1998).

6.2. A catalyst pellet problem We consider the problem of computing steady-state solutions to the nonlinear, reaction diffusion problem (c/#t = 92c− f 2c 2 in a long, cylindrical catalyst pellet. The catalyst pellet outer edge boundary condition is c(1)= 1 and we assume solutions are symmetric at the pellet center. Defining c(r)= u(r)+ 1, we obtain 0=

 

1 ( (u r − f 2(u+1)2 r (r (r

(14)

and boundary conditions #u(0)/#r = 0 and u(1)=0.

tfcts=gdf(x,’x.^(2*n)’,’x’,{’n’,[0:N]}); pts =colpts(tfcts,x,w.*(1-x.^2),[1],dx); [Q W A B] =colmat(tfcts,pts,x,w,dx,ddx); This example represents the case of slab geometry. w, A, B, and pts are the collocation-grid weight array, first-order differentiation array, the Laplacian operator

6.2.1. Collocation based on orthogonal polynomials We represent u(r, t) by the truncated trial function expansion u N(r)= cb. In this notation, the vector of trial functions and the corresponding mode amplitude coefficients are defined by c(r)= (1− r 2)[h1(r 2), …, hN (r 2)]

1056

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

b= [b1,b2, …, bN ]T.

c=tfcts*b+1;

By using only even powers of r, the trial functions cn satisfy the boundary condition at r = 0; the factor (1− r 2) forces the trial functions to also satisfy the boundary condition at r =1. The polynomials hn are constructed as the normalized Jacobi polynomials defined by a sequence of 2n-th order polynomials orthonormal with respect to weighted inner product Žhi, hj w = 10hihj (1− r 2)rdr =1 for i= j and Žhi, hj w = 0 for i" j, with h1 = 2. Choosing N=4 interior collocation points, we compute the collocation discretization array B and the discrete transformation array Q with colmat.m, after computing the collocation points with colpts.m. This allows us to write the discretized equations as

resid=ddr*c-(phi^2)*c.^2;

0= Bu− f 2(u+1)2.

(15)

The computational steps for generating the discretization arrays are [r w dr ddr] =pd(’disk’,ndisc); pr = 1-r.^2; f = gdf(r,’r.^(2*n)’,’r’,{’n’,[0:N]}); tfcts=f(:,1:N).*repmat(pr,1,N); clpts=colpts(f,r,w.*pr,[ ],dr); [Q W B] =colmat(tfcts,clpts,r,w,ddr); The solution at the collocation points is found from the values u that drive the collocation-discretized residual function to zero at the collocation points: colloc – resid=B*u-phi^2*(u + 1).^2; using a Newton–Raphson iteration technique. The residual function can be recovered from the discretized solution u by b = Q¯u;

Both the interpolated solution c and the residual function are plotted in Fig. 14. The convergence behavior as a function of interior collocation point number N is shown in Fig. 14 (right).

6.2.2. The stiff catalyst pellet problem A numerically stiff catalyst pellet problem (f= 1000) has been discussed in collocation literature (Villadsen & Michelsen, 1978; Sorensen & Stewart, 1982; Denison & Hamrin, 1983; Adomaitis & Lin, 1998) as a test case for evaluating different collocation procedures. Attempts at computing a solution in the slab geometry by a global polynomial collocation method were shown to fail (Villadsen & Michelsen, 1978) for discretizations performed using a moderate number of collocation points (e.g. N= 8). A spline collocation method, however, was shown to produce convergent results, measured in terms of the convergence of the effectiveness factor (Villadsen & Michelsen, 1978). Other BVP solution procedures have also been shown to be effective for solving this problem: Sorensen and Stewart (1982) demonstrated an exponential collocation method while Denison and Hamrin (1983) compared two BVP software packages based on finite differences and spline collocation. In this subsection, we demonstrate that this problem can, in fact, be solved with collocation techniques based on globally defined polynomial trial functions. The results are then compared to simulations performed using exponential trial functions (Guertin, Sorensen & Stewart, 1977; Sorensen & Stewart, 1982). The implementation procedure for these high-degree orthogonal polynomial methods in a slab geometry is identical to the methods discussed in the previous section except for the method of generating the orthogonal polynomials. When N is large, the high-degree functions x n cannot be evaluated accurately by gdf.m due

Fig. 12. Chemical process axial gas compression system illustrating the inlet and outlet duct and compression sections (left) and the compressor characteristics (right).

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1057

Fig. 13. Simulation results showing the growth of an infinitesimal spatial disturbance to a fully developed rotating stall cell in mode amplitude space (left). Snapshots of the rotating stall cell (right) illustrating the rotation of the traveling wave solution. Solutions correspond to N =21.

to the ill-conditioned Vandermonde matrix. Chebyshev polynomial formulations, however, can be used efficiently to generate the basis function h(x). The shifted Chebyshev formula can be used with x 2 substituted for x, i.e. Tn (x)=cos(n cos − 1(2x 2 −1)). The MATLAB command line to generate the trial functions is f = gdf(x.^2,’cos(n.*acos2*x-1))’ … ’x’,{’n’,[0:N]}); The functions are then orthonormalized with respect to (1−x 2) using gs.m. The results for the boundary flux (x= 1) as a function of polynomial truncation number is shown in Fig. 15 (case OC). We conclude that solutions can be computed using a global collocation technique, although the procedure becomes computationally intensive for large N. The problem also can be solved directly with a pseudospectral method, using the fine grid of the MWRtools. In other words, xˆ and B. computed from pd.m are used for the collocation points and discretization array as part of the orthogonal collocation method. However, the mixed collocation method must be employed and residual computations cannot be performed accurately in this single-grid mode. Representative results generated by this procedure are presented in Fig. 15 as the case PS showing the convergence of the boundary flux; a value of 816.4951 is found for N= 200, which corresponds to a 1.81×10 − 4% error. The difference between the OC and PS solutions in the range of N where the results overlap is attributable to the interior and mixed collocation formulations, respectively of the two methods. Finally, a collocation procedure based on exponential trial functions (Guertin et al., 1977; Sorensen & Stewart, 1982) was used that provides rapid convergence of the boundary flux. The procedure for generating the trial functions and collocation point locations can be reproduced using MWRtools using the following procedure.

The trial functions suggested in Sorensen and Stewart (1982) are (1− x) j − 1 (x − 1)f (1+ x) j − 1 − (x + 1)f e + e cj = 2 2 for j] 1, which can be generated with gdf.m. A trial function expansion of the form N

c(x)=2c1 + % aici (x) j=2

is used in conjunction with an interior collocation method. The location of collocation points is based on the orthogonality condition:

&

1

Q(x)ck (x)dx = 0

k= 1, …, N−1

0

and Q(x) is defined as N−1

Q(x)= cN (x)+ % cj (x)bj. j=1

The MWRtools commands for computing the bj and collocation points are b=-wip(psi(:,1:N-1),psi(:,1:N-1),w) … ¯wip(psi(:,N),psi(:,1:N-1),w); Q=psi(:,N) +psi(:,1:N-1)*b; clpts=nonzeros(rdf(x,Q,dx*Q)); We use the nonzeros on the output produced by rdf.m because the amplitude of Q(x) is vanishingly small and the zero roots are spurious (alternatively we could specify a smaller error stopping criterion in rdf.m). The solution based on a six-trial function expansion is shown in Fig. 16; the corresponding dimensionless boundary flux is 816.5823 (−1.05× 10 − 2% error). The solution and the boundary flux value agree with values reported in Sorensen and Stewart (1982). However, it is interesting to note that the residual between collocation points is still substantial (Fig. 16, right).

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1058

Fig. 14. Steady-state solution and the residual function (left) produced by the orthogonal collocation discretization technique for the catalyst pellet problem with f =4 and N = 4. The convergence behavior in terms of number of interior collocation is also shown (right).

d2f d2c and dc= dx 2 dy 2

7. Higher-dimensional projections

lf=

To introduce applications of the MWRtools library to problems with two spatial dimensions, consider the horizontal, thin, steel cooling fin shown in Fig. 17. The two-dimensional modeling equation for the dimensionless temperature u(x, y, t) developed under the thinvane assumption is

subject to the homogeneous boundary conditions (Eq. (17)) and u= 0 at y= 0. Computing these eigenfunctions and eigenvalues using

( 2u (u ( 2u = 2 +a 2r 2 −hu (y (t (x

[lambda,phi] =sl(’slab’,dx,x,1, …

(u =bu (x

n   x=0

[u = 1]y = 0

(u = −bu (x

(u = − bu (y

[y,wy,dy] =pd(’slab’,nd – y);

(16)

and boundary conditions



[x,wx,dx] =pd(’slab’,nd – x);

n

-beta,1,beta,wx); [delta,psi] =sl(’slab’,dy,y,0, …

n

1,1,beta,wy); (17)

x=1

y=1

with h= 0.2 and b =0.9 corresponding to physically relevant dimensionless parameter values. The aspect ratio is given by ar =Lx /Ly =0.5 for this problem. The dimensionless temperature is defined by u= (T −Ta )/ (Tr − Ta ), where Tr and Ta correspond to the root and ambient dimensional temperatures, respectively. The initial condition is given by u(x, y, 0) = 1; therefore this problem represents a sudden change in ambient temperature from a value equal to the root temperature. The solution to Eq. (16) subject to Eq. (17) and initial condition is represented by the combination of the two sets of trial functions

results in sets of orthonormal eigenfunctions. The orthogonality of the eigenfunctions is defined by the inner products  f(y),g(y)(y) =

&

1

f(y)g(y)dy

0

and an analogous inner product for the x direction. Substituting Eq. (18) into Eq. (16) and projecting this residual function onto fpcq gives

u(x,y,t)= uV(x,y,t)+u(V(x,y) N

N

i, j = 1

i=1

= % ai, j (t)fi (x)cj (y) + % bifi (x)(1 −y)2. (18) The trial function components, fi (x) and cj (y), are computed as nontrivial solutions to the eigenvalue problems

Fig. 15. Boundary flux in terms of number of trial functions used for the steady state catalyst pellet problem for the case of Thiele modulus f =1000 with slab geometry, comparing orthogonal polynomial collocation (OC) and pseudospectral discretization results (PS).

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

1059

Fig. 16. Steady state solution (left) for the catalyst pellet problem using exponential trial functions for the case of f=1000 and N= 5 collocation points. The plot of the residual function (right) with collocation positions marked.

dap,q =(lp + a 2r dq −h)ap,q +fp,q dt

for q=1:Q u=u+a(p,q)*phi(:,p)*psi(:,q)’;

= gp,qap,q +fp,q

end

by defining fp,q =[(1 − y)2(lp − h) + 2a 2r ]bp,cq (y)

end

and bn =Ž1,fn (x). Therefore,

Representative results are presented in Fig. 18, where the initial, final, and two snapshots of the transient temperature profile are shown.



ap,q (t)= ap,q (0)+

n

fp,q gp,q t fp,q e − gp,q gp,q

with

8. Concluding remarks

ap,q (0)= bp 2y −y 2,cq (y). Having computed the discretized eigenfunctions and eigenvalues, the coefficients resulting from projections of the initial conditions and inhomogeneous terms are computed using the weighted inner product function wip.m: b = wip(ones(nd – x,1),phi,wx); f = wip(psi,((1-y).^2*)(lambda’-h) … +2*alpha – r^2).*repmat(b’,nd – y,1),wy); a0 = b*wip(2*y-y.^2,psi,wy)’;

In this paper, we presented a set of computational tools for performing elemental steps in weighted residual method solution procedures and demonstrated the effectiveness of the MATLAB-based algorithms by solving a range of chemical engineering modeling problems. The development of a common set of computational techniques applicable to eigenfunction expansion, Galerkin, and collocation solution procedures was made possible by the discretized function representation structure employed in our computational framework. This approach also made possible accurate, and sometimes exact, assessment of the discretization error with

P = length(lambda); Q = length(delta); gamma=repmat(lambda,1,Q) + … alpha – r^2*repmat(delta’,P,1)-h; The mode amplitude coefficients ap, q are then computed and used to reconstruct the solution u(x, y, t) at selected points in time: a =(a0+f./gamma).*exp(gamma*time) … -f./gamma; u =(phi*b)*((1-y).^2)’; for p=1:P

Fig. 17. The geometry of the cooling fin. Temperature profiles for the fin u(x, y, t) in the x – y plane are to be computed with an eigenfunction expansion method.

1060

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

Fig. 18. Snapshots in time of the cooling fin temperature profile computed with P= Q = 20 eigenfunctions.

little additional computational effort. The functions and example scripts can be downloaded from the MWRtools website http://www.ench.umd.edu/ software/MWRtools. It is this combination of flexibility and the ability to quantify solution accuracy that has made it possible to use these numerical techniques as a powerful tool for developing and implementing advanced MWR, since advanced MWR applications can require a combination of the eigenfunction expansion, Galerkin, and collocation methods. For example, in Chang and Adomaitis (1999) a Galerkin discretization method was used to compute the gas flow field of a commercial chemical vapor deposition reactor which was subsequently used in a combined collocation/eigenfunction expansion solution of the gas temperature field to determine the wafer/gas phase heat transfer coefficient. In Adomaitis and Lin (1999), an exact discrete analog to the Galerkin projection was developed to study the convergence behavior of the discretized heterogeneous catalysis problem in a cylindrical geometry (discussed in Section 6.2.1) as a function of Thiele modulus. The authors found that the convergence problems of loworder discretized solutions were attributable a saddlenode bifurcation between the true solution and a

spurious solution. An application of a high-degree, pseudospectral discretization procedure based on globally-defined trial function expansions was presented in Lin and Adomaitis (1998) for simulation of a DC glow discharge system. Model reduction methods can also be developed in the MWRtools framework. In Chang and Adomaitis, (1998a,b), a collocation equivalent to a nonlinear Galerkin procedure was presented as a method for reducing the number of dynamic degrees of freedom in a reduced-order chemical vapor deposition system simulation. A different model reduction approach, based on identifying the dominant wafer temperature modes and using these optimized trial functions as part of a collocation discretization, was developed in Theodoropoulou et al. (1998) to produce a reduced-order model of a rapid thermal processing system. Enhancements to the MWRtools library currently under development include efficient algorithms for projections in two- and three-dimensional physical domains, collocation-based ODE-algebraic equation system numerical integrators, and other computational method improvements. Updated information will be made available at the MWRtools website.

Y.-h. Lin et al. / Computers and Chemical Engineering 23 (1999) 1041–1061

Acknowledgements This work was supported by the National Science Foundation through grant NSF EEC 95 – 27 576 and the Petroleum Research Fund of the American Chemical Society through grant 31391-G9. The authors would like to thank Dr H. S. Brown for his constructive comments on an early draft of this manuscript and the contributions of Dr A. Theodoropoulou, Dr A. J. Newman, and Mr P. Loezos to the preliminary version of the MWRtools library.

References Adomaitis, R. A., & Abed, E. H. (1996). Bifurcation analysis of nonuniform flow patterns in axial-flow gas compress. In V. Lakshmikantham, Proceedings of the first world congress of nonlinear analysts 1992 (pp. 1597–1608). Adomaitis, R. A., & Lin, Y.-h. (1998). A technique for accurate collocation residual calculations. Chemical Engineering Journal, 71, 127 – 134. (Also in ISR TR 98–6.) Adomaitis, R. A., & Lin, Y.-h. (1999). A collocation/quadraturebased Sturm – Liouville solver. Applied Mathematics and Computation. (Accepted for publication; also in ISR TR 99–1). Articolo, G. A. (1998). Partial differential equations and boundary 6alue problems with Maple V. San Diego, CA: Academic. Ascher, U., Christiansen, J., & Russell, R. D. (1981). Collocation software for boundary-value ODEs. ACM Transactions on Mathematical Software, 7 (2), 209–222. Carberry, J. J. (1976). Chemical and catalytic reaction engineering. Chemical engineering series (pp. 502–503). New York: McGraw Hill. Chang, H.-Y., & Adomaitis, R. A. (1998a). Model reduction for a tungsten chemical vapor deposition system. In Proceedings of DYCOPS-5 (pp. 35 –40). Corfu, Greece. Chang, H.-Y., & Adomaitis, R. A. (1998b). Model reduction for tungsten chemical vapor deposition. In Proceedings of the 1998 American Control Conference (pp. 1365–1366). Philadelphia, PA. Chang, H.-Y., & Adomaitis, R. A. (1999). Analysis of heat transfer in a chemical vapor deposition reactor: an eigenfunction expansion solution approach. International Journal of Heat and Fluid Flow, 20, 74 – 83. Childs, B., Scott, M., Daniel, J. W., Denman, E., & Nelson, P. (1979). Codes for boundary-6alue problems in ordinary differential equations. Lecture notes in computer science, no. 76. New York: Springer-Verlag. Cooper, J. M. (1998). Introduction to partial differential equations with MATLAB. Boston, MA: Birkhauser. Denison, K. S., & Hamrin, C. E. Jr. (1983). Solution of boundary value problems using software packages: DD04AD and COLSYS. Chemical Engineering Communications, 22, 1–9. Don, W. S., & Solomonoff, A. (1997). PseudoPack-pseudo-spectral differentiation software package user manual for 6ersion 2.3 beta. Providence, RI: Brown University. Eaton, W. J. (1995). Octa6e—a high-le6el interacti6e language for numerical computations, edition 1.1 for Octave version 1.1.1. Finlayson, B. A. (1972). The method of weighted residuals and 6ariational principles. New York: Academic.

.

1061

Finlayson, B. A. (1980). Nonlinear analysis in chemical engineering. New York: McGraw-Hill. Gay, D. H., & Ray, W. H. (1995). Identification and control of distributed parameter systems by means of the singular value decomposition. Chemical Engineering Science, 50, 1519 –1539. Gottlieb, D., & Orszag, S. A. (1977). Numerical analysis of spectral methods: theory and applications. No. 26 Philadelphia, PA: Society for Industrial and Applied Mathematics. Gottlieb, D., & Shu, C.-W. (1997). The Gibbs phenomenon and its resolution. SIAM Re6iew, 39, 644 – 668. Guertin, E. W., Sorensen, J. P., & Stewart, W. E. (1977). Exponential collocation of stiff reactor models. Computers and Chemical Engineering, 1, 197 – 203. Lin, Y.-h., & Adomaitis, R. A. (1998). A global basis function approach to DC glow discharge simulation. Physics Letters A, 243 (3), 142 – 150. (Also in ISR TR 97 – 81.) MacCluer, C. R. (1994). Boundary 6alue problems and orthogonal expansions —physical problems from a Sobole6 6iewpoint. New York: IEEE. Madsen, N. K., & Sincovec, R. F. (1979). Algorithm 540: PDECOL, general collocation software for partial differential equations. ACM Transactions on Mathematical Software, 5, 326 – 351. MathWorks (1995). Partial differential equation toolbox-for use with MATLAB User’s Guide. Natick, MA. Michelsen, M. L., & Villadsen, J. (1972). A convenient computational procedure for collocation constants. Chemical Engineering Journal, 4, 64 – 68. Michelsen, M. L., & Villadsen, J. (1981). Polynomial solution of differential equations. In Mah, R. S. H., & Seider, W. D., Foundations of computer-aided chemical process design, Proceedings of an international conference (pp. 341 – 368). New York: Engineering Foundation. Moore, F. K., & Greitzer, E. M. (1986). A theory of post-stall transients in axial compression systems: part I-development of equations. ASME Journal of Engineering for Gas Turbines and Power, 108, 68 – 76. Ogunnaike, B. A., & Ray, W. H. (1994). Process dynamics, modeling, and control. New York: Oxford University. Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. (1988). Numerical recipes in C. New York: Cambridge University. Ray, W. H. (1981). Ad6anced process control. Chemical engineering series. New York: McGraw-Hill. Rice, R. G., & Do, D. D. (1995). Applied mathematics and modeling for chemical engineers. New York: J. Wiley & Sons. Sorensen, J. P., & Stewart, W. E. (1982). Collocation analysis of multicomponent diffusion and reactions in porous catalysts. Technical report, University of Wisconsin Mathematics Research Center, Madison, WI. Stroud, A. H. (1974). Numerical quadrature and solution of ordinary differential equations. In Applied mathematical science series, vol. 10. New York: Springer-Verlag. Theodoropoulou, A., Adomaitis, R. A., & Zafiriou, E. (1998). Model reduction for optimization of rapid thermal chemical vapor deposition systems. IEEE Transactions on Semiconductor Manufacture, 11 (1), 85 – 98. Villadsen, J., & Michelsen, M. L. (1978). Solution of differential equation models by polynomial approximation. International series in physical and chemical engineering science. Englewood Cliffs, NJ: Prentice-Hall Villadsen, J. V., & Stewart, W. E. (1967). Solution of boundary-value problems by orthogonal collocation. Chemical Engineering Science, 22, 1483 – 1501.