Available online at www.sciencedirect.com

Computers and Chemical Engineering 32 (2008) 1305–1315

On POD reduced models of tubular reactor with periodic regimes Katarzyna Bizon a,∗ , Gaetano Continillo a , Lucia Russo b , Joanna Smuła c a

b

Dipartimento di Ingegneria, Universit`a del Sannio, Piazza Roma 21, 82100 Benevento, Italy Dipartimento di Ingegneria Chimica, Universit`a Federico II, Piazzale Tecchio, 80, 80125 Naples, Italy c Wydział Chemiczny, Politechnika Sl˛ ´ aska, ul. Strzody 9, 44-100 Gliwice, Poland Received 27 January 2007; received in revised form 2 June 2007; accepted 6 June 2007 Available online 15 June 2007

Abstract The distributed model of a tubular reactor with recycle is reduced by approximating it first with a Continuous Stirred Tank Reactor (CSTR) cascade. Proper orthogonal decomposition (POD) with Galerkin projection is introduced to study oscillatory regimes. The dynamics of the resulting models is studied via numerical simulation, and solutions in the form of time series and phase plot diagrams are compared to those obtained for the original CSTR cascade model. Different methods to choose the minimum number of basis functions are compared and discussed. Solution diagrams built with POD models are compared with those from CSTR cascade, as a function of the Damk¨ohler number. Features and limitations of POD models are discussed for different snapshot sampling policies and for different values of the P´eclet number. Qualitative performance increases and quantitative performance decreases as samples from steady states are considered along with periodic solution samples. Performance becomes worse as the P´eclet number increases. © 2007 Elsevier Ltd. All rights reserved. Keywords: Model reduction; Proper orthogonal decomposition; Non-linear analysis; Tubular chemical reactor with recycle; Bifurcation analysis

1. Introduction Many chemically reactive systems are characterized by parametric sensitivity that results in the occurrence of a wide variety of static and dynamic phenomena (Aris, Aronson, & Swinney, 1991). Depending on the operating conditions or on the values assumed by different parameters of the model which describes the physical and chemical properties of the system, static equilibrium as well as more complex time-asymptotic regimes such as periodic, quasi-periodic or chaotic oscillations can be observed. To describe all possible states and determine values of parameters for which dramatic qualitative or quantitative changes of the solution occur, it is necessary to conduct an accurate bifurcation analysis in the parameter space. In order to accomplish this, the original PDE evolutionary equations, which in most cases constitute the model, are to be reduced to a finite-dimensional dynamical system. This is typically done ∗

Corresponding author. Present address: Istituto di Ricerche sulla Combustione CNR, Naples, Italy. Tel.: +39 081 5704105; fax: +39 081 5936936. E-mail addresses: [email protected] (K. Bizon), [email protected] (G. Continillo), [email protected] (L. Russo), [email protected] (J. Smuła). 0098-1354/$ – see front matter © 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2007.06.004

by numerical discretization of the spatial differential operators with finite difference schemes or various spectral methods (see, for example, Quarteroni & Valli, 1994, chap. 5). The resulting reduced model can then be analyzed in its dynamical behaviour via bifurcation analysis (Kuznetsov, 1998). In spite of the availability of various packages for performing bifurcation analysis (Doedel et al., 2001), often the problem becomes very difficult because the order of the dynamical model required to properly describe the system is very high. This, coupled with obvious limitations of software and hardware, leads to the need of expressing the original model as a set of ordinary differential equations of the lowest possible order. Among classical numerical approaches, finite difference methods are simple but require a relatively large number of ODE as compared, for instance, to orthogonal collocation (Villadsen & Michelsen, 1978). Other spectral methods have been used for reactive systems such as Galerkin projection (see, for example, Oran & Boris, 1987). Projection methods provide an interesting framework in that, by proper choice of the functional basis, one can reduce the number of ODEs necessary to accurately describe the dynamics of the original PDE model. Many applications adopt as functional basis a set of eigenfunctions of the main differential operator in the model equations (Adrover, Cerbelli, & Giona, 2002) which by

1306

K. Bizon et al. / Computers and Chemical Engineering 32 (2008) 1305–1315

Nomenclature Da f Le n PeH PeM

Damk¨ohler number recycle coefficient Lewis number order of reaction P´eclet number for heat P´eclet number for mass

Greek letters α conversion degree β coefficient related to adiabatic temperature increase δ heat transfer coefficient γ activation energy θ dimensionless temperature τ dimensionless time ξ dimensionless position coordinate along reactor Subscripts and superscripts out reactor outlet α conversion degree θ dimensionless temperature 0 refers to initial profile

the way in our case (Laplacian expressing molecular diffusion) is a self-adjoint operator and thus it has a set of orthonormal eigenfunctions. The number of basis functions chosen for the approximation is rather critical: Graham and Kevrekidis (1996) showed that an insufficient number of basis functions – in their case Chebyshev polynomials – can cause that significant phenomena such as period doubling are not detected. Hence, they proposed to use basis determined by proper orthogonal decomposition, with some extension of the original procedure that yield optimal basis for parameter varying applications. Proper orthogonal decomposition (POD) is a procedure that delivers an optimal set of empirical basis functions from an ensemble of observations obtained either experimentally or from numerical simulation, which characterize the spatio-temporal complexity of the system (see, for example, Holmes & Lumley, 1996). Obtained orthogonal functions can be afterwards used in a Galerkin projection of the original system and as a result a low-dimensional model can be developed (Lumley, 1967). POD has been widely used in disciplines such as control of chemical reactors (Padhi & Balakrishnan, 2003; Shvartsman et al., 2000), modelling of nonreactive fluidized beds (Cizmas, Palacios, O’Brien, & Syamlal, 2003; Yuan, Cizmas, & O’Brien, 2005), as well as in reducing complex chemical kinetic models (Danby & Echekki, 2005). Besides, it has been applied in many engineering fields, among other things to damage detection (Joyner, 2004), low-order modelling of flows (Sirsup, Karniadakis, Xiu, & Kevrekidis, 2005) or microelectromechanical systems (Liang, 2002) and, depending on the field of application, it is also known as Principal Components Analysis (PCA) (Jolliffe, 2002), Sin-

gular Value Decomposition (SVD) (Grossman, Kamath, & Kegelmeyer, 2001, chap. 22), Karhunen–Lo`eve (KL) decomposition (Kirby & Sirovich, 1990) or Hotelling transformation (Eriksson, Jim´enez, B¨uhler, & Murtagh, 2002). Some extensions of the procedure, such as application of POD to moving boundary problems of fluid flow was studied by Utturkar, Zhang, and Shyy (2004), whereas Gokulakrishnan, Lawrence, McLellan, and Grandmaison (2006), introduced a functional approach (functional-PCA) with application to oxidation processes. In this paper, we consider a model of pseudohomogeneous tubular reactor with mass recycle. This system exhibits the socalled “relaxation oscillations” (Phillipson & Schuster, 2001; Smuła, 2006), an oscillatory regime in which a fast variation of the state is followed by a slow “relaxation” phase. Such type of regime is difficult to handle from a numerical point of view. Moreover, periodic regimes retain both fast and slow modes. These two features identify this system as an appropriate benchmark for testing reduced models. We apply POD/Galerkin and analyze the performance of the method by comparing solutions from the reduced model with a “reference” numerical solution. 2. Mathematical model of the reactor The model of pseudohomogeneous tubular reactor with mass recycle is essentially that described in Berezowski (2000). The dimensionless mathematical model is given by the following system of mass and heat balance equations:   ∂α ∂α βθ 1 ∂2 α n + (1 − f )Da(1 − α) exp γ + = ∂τ ∂ξ PeM ∂ξ 2 1 + βθ (1) Le

∂θ ∂θ 1 ∂2 θ + (1 − f )Da(1 − α)n + = ∂τ ∂ξ PeH ∂ξ 2   βθ × exp γ + (1 − f )δ(θH − θ) 1 + βθ

with the boundary conditions resulting from recycle:  1 ∂α  ; fα(1) = α(0+ ) − PeM ∂ξ ξ=0+  1 ∂θ  + fθ(1) = θ(0 ) − PeH ∂ξ ξ=0+  ∂α  = 0; ∂ξ ξ=1

 ∂θ  =0 ∂ξ ξ=1

(2)

(3)

(4)

where α is a degree of conversion, θ the dimensionless temperature and 0+ denotes right-side of inlet of the reactor. This system was also studied earlier by Liu and Jacobsen (2004) with a reduced model based on orthogonal collocation incorporating error estimates and dynamic mesh control. In the present work, the PDE model is first reduced to a system of ODE by approximating it with a cascade of N CSTRs (Levenspiel, 1998). The Continuous Stirred Tank Reactor (CSTR) cascade is then employed to build reference solutions

K. Bizon et al. / Computers and Chemical Engineering 32 (2008) 1305–1315

by time integration, to conduct bifurcation analysis, and finally to collect data needed for the generation of POD basis functions with a procedure described below. The number N of tanks in the cascade is chosen equal to 1/2Pe, which can be shown to be an optimal approximation basing on Residence Time Distribution (RTD) arguments, as shown, for example, in Bretsznajder, Kawecki, and Marcinkowski (1973). 3. Proper orthogonal decomposition

1307

Using the POD modes, the solution ut (x) can be expressed by a linear combination of the eigenfunctions: u˜ t (x) =

K

ak (t)ϕk (x)

(9)

k=1

where K < N is the number of modes used for truncation, whereas ak (t) are modal coefficients that can be calculated by Galerkin projection of the original PDE on the POD modes. Namely, consider a PDE of the form (Rowley, Colonius, & Murray, 2001):

3.1. Theory

u˙ = Dμ (u)

In the POD scheme, the objective is to determine a set of orthogonal basis functions which minimize, on average, the least square error between the truncated representation of the model and the true solution. Suppose we have given a time series, obtained from simulation or experiment, ut (x), where t denotes time and x denotes position in space. In practice, in numerical models the spatial domain is discretized and the number of time samples is finite, therefore usually the sampled data set is a vector-valued function given as a matrix: ⎡ ⎤ u1 (x1 ) u2 (x1 ) · · · uM (x1 ) ⎢ u (x ) u (x ) · · · u (x ) ⎥ 2 2 M 2 ⎥ ⎢ 1 2 ⎥ U=⎢ (5) .. .. ⎢ .. ⎥ .. . ⎣ . ⎦ . .

where u = ut (x) and Dμ (·) is a non-linear operator that involves spatial derivatives and depends on some set of parameters μ. Replacing u by a truncated linear combination of the first K POD modes (10) and conducting projection onto these modes, we obtain the following set of ODEs written in terms of the modal coefficients:

u1 (xN )

u2 (xN )

· · · uM (xN )

where N is the number of positions in the spatial domain and M is the number of samples taken in time. It can be shown that a suitable POD basis Φ = {ϕ1 , ϕ2 , . . ., ϕN } is obtained by solving the eigenvalue problem: CΦ = λΦ

(6)

where C is the averaged autocorrelation matrix: C(x, x ) = Ut (x), Ut (x )

(7)

and angular brackets denote time-averaging operation. Besides, it is common in POD applications to subtract the time average from each column of the matrix of initial data U, although, as mentioned in Crommelin and Majda (2004), this operation does not improve the performance of POD models: simply, if the time average is not subtracted when calculating C, the leading POD mode coincides with the time-average state. The ordering of the eigenvalues from the largest to the smallest induces an ordering in the corresponding eigenfunctions, from the most to the least important. Hence, in order to determine the truncation degree of the POD reduced model, we define the cumulative correlation energy captured by the K successive modes which is given by: K λi EK = k=1 . (8) N k=1 λi Thus, if we set the minimum required cumulative correlation energy the number of necessary modes is determined according to Eq. (8).

a˙ k = (Dμ (˜u), ϕk ),

(10)

k = 1, . . . , K

(11)

where (·, ·) denotes inner product, while initial values can be obtained by projecting u0 (x) on to eigenfunctions. The coefficients ak (t) can be also computed using linear (Joyner, 2004) or cubic spline interpolation (Bui-Thanh, Damodaran, & Willcox, 2003), based on the known values of the coefficients calculated earlier for instants of time or values of parameter in which samples were taken. It is important to notice that such procedure does not provide a reduced system of ODEs that could be for example used in continuation method or system control; it just provides a means for the calculation of solutions for intermediate values in the range of conditions (such as time of integration, range of parameters) for which the original, full order system was solved. Moreover, for the reliability of the interpolation, smoothness of the coefficients is required: this means that interpolation may not function when used, for example, across a catastrophic bifurcation point. Nevertheless, POD/Interpolation can be used for very fast computation of the coefficients, starting from the values of the coefficients known from calculations conducted for different parameter values. The method allows reconstructing continuous space–time (or parameter) description of the solution. Hence, it can be successfully applied to reconstruct the behaviour of a system from experimental data, when a mathematical model is not available (Druault, Guibert, & Alizon, 2005). 3.2. Sampling The reliability of the POD basis, i.e., the ability of reduced models to capture correctly the dynamics of the system, depends on the ensemble of observation, that is, in our case, on the generation of this ensemble. Time of sampling, location and number of samples may be arbitrarily chosen in view of constructing the best possible basis. Moreover, if the system model depends on one or more parameters, the reliability of the basis as the parameters change will depend on how values of the parameters are assigned in generating the samples. Finally, initial conditions can be changed either to best capture transient behaviour or, in

1308

K. Bizon et al. / Computers and Chemical Engineering 32 (2008) 1305–1315

case of multiple solutions, to characterize all different regimes (Alonso, Kevrekidis, Banga, & Frouzakis, 2004). In a recent study of dynamical models for control of cavity oscillations, two sets of POD modes were built (Rowley et al., 2001) from samples taken in the transient as well as in the fully developed region. The behaviour of POD/Galerkin models for different initial conditions was compared. It was reported that modes constructed from samples taken from the oscillatory dynamic regimes work much better for long time behaviour than those built from samples taken from the transient. Moreover, the authors showed that including samples from the transient increases the number of modes required for the determination of an accurate reduced order model or can even produce instabilities where the full model is stable. It was also mentioned that, as one can expect, such a model is not valid far from the parameter range in which the snapshots were taken. Adrover and Giona (2003) address a sampling issue in illustrating an alternative approach to classical POD, in which the PDE is projected onto a linear subspace spanned by so-called Snapshot Archetypes. These latter are selected, among available solutions of the PDE, by maximizing some mutual distance, in view of optimizing the span of this subspace over the solution set for the original PDE. Non-linear systems may obviously exhibit multiple solutions for a single chosen set of parameter values. These solutions may have different dynamical character (e.g. stationary steady states coexisting with oscillatory solutions) (Subramanian & Balakotaiah, 1997). The system trajectory towards one or the other depends on initial conditions. In order to predict system behaviour for a wider range of conditions, Graham and Kevrekidis (1996) propose some modification of the procedure for the construction of the data ensemble. Particularly, with the purpose of determining a POD basis that contains global information about the model, i.e., represents all its possible timeasymptotic states, the data ensemble was composed of 1000 initial conditions integrated forward in time. All computations were done for one given value of the parameter, for which chaos was known to exist. Interestingly, for this system, even though all these conditions were integrated at one value of the parameter, the obtained POD basis was able to approximate well the solutions manifold. On the other hand, it was shown (Zhang, Henson, & Kevrekidis, 2004) that, in order to capture the dynamics of a

system representing a cell population model over a wide range of conditions, it is necessary to incorporate in the data ensemble samples for several parameter values that span the desired region. Also, it was stated that the combined data set should contain samples from the transient as well as from the fully developed region, i.e., the final attractor. In the mentioned work, the POD basis describing the global behaviour of the system was built from data sets collected from simulations of the original system at different initial conditions and up to six different values of the bifurcation parameter. Then, the reduced model was successfully applied to conduct bifurcation analysis. 4. Analysis, results and discussion Our goal is to evaluate the performance of POD/Galerkin. To this aim, solutions obtained by this approach are compared with those obtained with the CSTR cascade model. Parameters of the model were kept constant as follows: f = 0.2, β = 0.75, Le = 5, γ = 15, δ = 0 and n = 1. Fig. 1a reports the solution diagrams for the system under study, obtained by parameter continuation of a 50 tank series approximation, for PeM = PeH = 100 and the Damk¨ohler number Da chosen as the bifurcation parameter. As it is seen, for the assumed above set of parameters, multiple solutions exist, i.e., one stable fixed point and stable and unstable oscillations. Here, recalling the above mentioned observations by Rowley et al. (2001), and the criteria stated by Adrover and Giona (2003), for optimal selection of their Snapshot Archetypes, we note that samples from periodic regimes are expected to better span the realm of possible space profiles. Therefore, we generate snapshots only for stable oscillating solutions. To this aim we chose suitable initial conditions and range of variation of the Damk¨ohler number Da ∈ [0.03, 0.08]. In order to determine the initial profiles in the bed of the reactor that for the whole assumed range of Da gives as a result oscillating behaviour, it would be necessary to compute the separation line between basins of attraction of the fixed point and the stable oscillations. This is a formidable problem when faced in the whole space state for a distributed system like ours. In practice, tubular reactors may be thought as systems starting from spatially uniform initial conditions and thus, with the purpose of simplification, we considered only classes of constant pro-

Fig. 1. (a) Solution diagram for conversion degree vs. Damk¨ohler number, CSTR cascade approximation. Stable stationary solutions and stable oscillations (—); unstable oscillations (– – –); unstable stationary solutions (· · ·). (b) Boundaries between attraction basins for different values of Damk¨ohler number.

K. Bizon et al. / Computers and Chemical Engineering 32 (2008) 1305–1315

1309

files of conversion degree α0 and dimensionless temperature, θ 0 . This assumption allows to calculate boundaries of the socalled reduced basins of attraction (Russo, Mancusi, Continillo, & Crescitelli, 2001), which in this case is a curve in the α0 , θ 0 space. The algorithm is based on the simple bisection procedure for root-finding. Starting from two initial conditions leading to different type of solutions, i.e., belonging to different basins of attraction of the system, the first point of the separation curve is determined by means of a halving algorithm along the line connecting these initial conditions. Then, successive points are determined similarly, operating on arcs rather than on intervals, in the circle of diameter ρ centred in the previous point of separatrix. 4.1. POD/Galerkin As we can see in Fig. 1b, if the initial temperature is sufficiently high, then for all values of Da in the assumed interval all solutions approach the stable limit cycle. Therefore, we chose the initial profiles as α(ξ, 0) = α0 = 0.2 and θ(ξ, 0) = θ 0 = 0.3 and integrate the system in time for Da = 0.03, 0.04, 0.05, 0.06, 0.07 and 0.08. Rowley et al. (2001) showed that developing modes, i.e., modes constructed from samples taken in the transient, can give worse result and make the model unstable but, on the other hand, Zhang et al. (2004) reported that it is necessary to include in the ensemble transient data as well as stationary data. Therefore, for each of value of Da we took 300 equally spaced temporal samples between τ = 0 and 150, that is starting from the initial condition until the trajectory of the solution settle down onto stable oscillations. Fig. 2 shows a typical sequence of snapshots of dimensionless temperature space profiles taken from the transient Da = 0.03. Following the POD procedure described above, collected data sets were combined into one matrix of snapshots for each state variable (Padhi & Balakrishnan, 2003) and manipulated to produce two sets of orthogonal basis functions ϕk (ξ), ψk (ξ), k = 1, . . ., N. This makes it possible to determine the order of approximation separately for each group of state variables. In fact, generally, a different number of modes is applied for each variable. Hence, state variables α(ξ, τ) and θ(ξ, τ) can be expressed as a linear combination of the basis functions: α

˜ τ) = α(ξ,

K k=1

θ

ak (τ)ϕk (ξ),

˜ τ) = θ(ξ,

K

bk (τ)ψk (ξ)

(12)

k=1

where Kα and Kθ are the – independently chosen – truncation ranks for conversion degree and dimensionless temperature, respectively, whereas ak (τ) and bk (τ) are time dependent coefficients of basis functions. 4.1.1. Correlation energy As a first application, ranks of truncation were chosen by calculating the corresponding values of the cumulative correlation energy until it reached a chosen value. Calculated amounts of energy and corresponding eigenvalues for leading modes, as well as some higher order modes that were used for model reduction, are collected in Table 1.

Fig. 2. Snapshots of temperature collected from the CSTR simulation, Da = 0.03 for dimensionless times τ from 0 to 9.

In order to build reliable reduced order models, the number of basis functions used in the projection should be chosen such as the cumulative energy of the modes EK be greater or equal to 99%. This criterion would lead to choose Kα = 4 and Kθ = 6 (Table 1). For such a value of cumulative correlation energy one might expect the model to be nearly perfect, however it was shown (Aubry, Lian, & Titi, 1993) that sometimes modes that represent even more than 99.999% might not reproduce the dynamics of the system satisfactorily. On the other hand, Joyner (2004) showed that taking too many modes can produce some noise-like behaviour or make the model unstable (Rowley et al., 2001). 4.1.2. Eigenvalue pairing Another problem appears when dealing, as it is in our case, with periodic behaviour. Rowley et al. (2001) showed that, whenever there are oscillations at a single frequency, eigenvalues occur in pairs. By “pairs” they intend two consecutive eigenvalues which differ in magnitude much less than the average

1310

K. Bizon et al. / Computers and Chemical Engineering 32 (2008) 1305–1315

Table 1 Largest eigenvalues of correlation matrices for conversion degree and temperature and their corresponding cumulative correlation energy No. of mode

1 2 3 4 5 6 7 10 11 12 13 18 19

Eigenvalues

Cumulative energy captured [%]

λαk

λθk

α EK

θ EK

23.2587 0.5683 0.2952 8.8597E−002 6.7855E−002 2.9643E−002 2.5002E−002 6.3135E−003 5.4795E−003 3.3460E−003 2.7619E−003 5.6317E−004 4.8125E−004

23.3869 2.7950 1.6347 0.3819 0.3031 0.1043 9.1178E−002 1.7742E−002 1.6118E−002 8.1701E−003 7.7015E−003 1.0386E−003 9.4812E−004

95.3882 97.7189 98.9299 99.2932 99.5715 99.6931 99.7956 99.9204 99.9429 99.9566 99.9680 99.9915 99.9935

81.0918 90.7833 96.4515 97.7758 98.8268 99.1888 99.5050 99.8338 99.8897 99.9180 99.9447 99.9881 99.9914

difference between two consecutive eigenvalues in the spectrum. This is especially appreciable when POD modes are determined from snapshots extracted from the long-term behaviour of the system. It was also found that low-order models can be unstable if both modes that correspond to each pair of eigenvalues are not kept: as a consequence, increasing the order of the model by incorporating only one additional eigenfunction could result into worse performance, even though we know that cumulative correlation energy increases, if the added eigenfunction belongs to a pair and the second eigenfunction of the pair is not incorporated as well. Hence, in this case the cumulative correlation energy appears to be not a sufficient criterion for the determination of the truncation order. Fig. 3a shows the full eigenvalue spectrum of POD modes for both state variables of our system model. Fig. 3b shows the decay of the first 10 eigenvalues. After the first dominant eigenvalues, pairing of eigenvalues 2, 3; 4, 5 and so on is still recognizable, even though the matrix does not correlate data collected from a simple oscillating system but rather incorporates data from simulations conducted for different values of the parameter, hence for different values of the oscillation period. We recall that, according to Perron’s theorem (Horn, 2002; Tarazaga, Raydan, & Hurman, 2001), any primitive matrix (i.e., a square non-negative matrix some power of which is positive) has a positive eigenvalue λ that is algebraically simple, equal

Fig. 3. (a) Eigenvalue spectrum for conversion degree () and temperature (). (b) Decay of the first 10 eigenvalues.

to the spectral radius and strictly greater than the modulus of any other eigenvalues. Moreover, this λ has associated a strictly positive eigenvector. Notice (Fig. 4a and b) that both the correlation matrix of conversion degree and the correlation matrix of temperature have strictly positive entries; hence, the presence of a dominant eigenvalue comes directly from the theorem.

Fig. 4. Autocorrelation matrices for conversion degree (a) and temperature (b).

K. Bizon et al. / Computers and Chemical Engineering 32 (2008) 1305–1315

Fig. 5. Leading POD basis functions for conversion degree (a) and temperature (b); first mode (—); second mode (– – –); third mode (· · ·).

Fig. 5a and b reports space diagrams of the three leading eigenfunctions. It appears that the first mode for each state variable is strictly positive, whereas the second and the third one attain negative values. In connection with this property, we observe that, often, an odd number of POD modes (the dominant one plus a number of pairs) gives a better reduced model than one containing one more eigenfunction. This is more visible for low-order modes. In order to verify the effectiveness of our global POD basis, we focused on performance of the low-order model for values of the Damk¨ohler number that were not included in the snapshots set. Here we report the case of Da = 0.035. In our case, projection of the original system onto four modes for conversion degree and six modes for dimensionless temperature (these modes represent 99% of the correlation energy of the system) should give satisfactory results. However, for these values the integrator even failed to converge, as a result of separation eigenfunctions that correspond to paired eigenvalues. Adding one higher mode for each state variable still yielded very inaccurate predictions of the solution. As we can see in Fig. 6, a reduced order model that captures a bit more than 99% energy of the system (5 + 7 modes) is unable to predict correctly either the transient or the long-term behaviour. Adding higher order modes improves significantly the solution. The model covering 99.9% of correlation energy (Fig. 7, 11 + 13 modes) yields results that are comparable with the full model in terms of amplitude and period of oscillations, however distinct disturbances in the relaxation interval of the oscillation of conversion degree cannot be accepted because the value of conversion degree significantly exceeds one. Fig. 8 shows results of the projection of the model onto 19 basis functions for each state variable. This number of modes represents 99.99% of the energy. As we can see on the plot, the 19-mode model makes the solution much smoother and also predictions of amplitude and period are more accurate. 4.1.3. Relative error versus energy Consideration of some norm of the error and not just energy seems to be essential in judging the quality of the approxima-

1311

Fig. 6. Comparison of time series for conversion degree (a) and temperature (b); POD/Galerkin with five modes for conversion degree and seven modes for dimensionless temperature (—) (99% of correlation energy captured); CSTR model (– – –).

tion. For instance, the relative error could be simply defined as the mean square error of the truncation ||u − u˜ ||2 , where ||·|| denotes Euclidean norm and · denotes time or space average. Fig. 9 shows the error of the POD approximation calculated for the interval of time from τ = 0 to 150, that contains transient as well as stable oscillatory regime, depending on the number of basis functions is shown. However, for steady oscillatory solutions, the presence of even small errors in the early transient could be misleading. In fact, errors in the transient can produce a phase shift, resulting in a large relative error. On the other hand, even though approximation in the transient is good and no phase shift is produced, the relative error can be affected by different problems, such as amplitude differences or frequency drift. Fig. 10 reports time series for 22-, 23- and 24-mode models. As we can notice on the plot of the error (Fig. 9), for these numbers of eigenfunctions a strange phenomenon occurs. Solu-

Fig. 7. Comparison of time series for conversion degree (a) and temperature (b); POD/Galerkin with 11 modes for conversion degree and 13 modes for dimensionless temperature (—) (99.9% of correlation energy captured); CSTR model (– – –).

1312

K. Bizon et al. / Computers and Chemical Engineering 32 (2008) 1305–1315

Fig. 8. Comparison of time series for conversion degree (a) and temperature (b); POD/Galerkin with 19 modes both for conversion degree and dimensionless temperature (—) (99.99% of correlation energy captured); CSTR model (– – –).

tion approximated by using 22 modes produce larger error than the solution of the 21-mode model. Adding one higher mode causes a distinct drop of the error (23-mode model), but then, the subsequent mode increases it back (24-mode model). As we can see in Fig. 10a, a slight error in the frequency produces an increasing error, up to one-half period shift, over times that are much larger than the minimum period (Fig. 10b). Next, in Fig. 10c and d we can see that the prediction of the frequency of the oscillations from the 23-mode model is much better; however, adding next mode (Fig. 10e and f) again produces a similar error as the one for the 21-mode model. This explains why the error defined earlier exhibits erratic behaviour as the number of modes is increased. 4.1.4. Approximation of the limit set Another important evaluation for the quality of the reduced model is the approximation of the limit set, namely how close to the limit set of the “reference” CSTR solution a reduced model can get. To evaluate the distance between limit sets we use the

Fig. 10. Comparison of time series in the transient (left columns) and at steady state (right columns) for temperature from the 22-mode (a and b), 23-mode (c and d) and 24-mode models (e and f); POD/Galerkin (—); CSTR model (– – –).

Hausdorff distance between sets defined as follows (Stuart & Humphries, 1998): dH (A, B) = max{dist(A, B), dist(B, A)}

(13)

where dist(B, A) = sup dist(u, A) u∈B

(14)

is called the Hausdorff semi-distance, and: dist(u, A) = inf u − v v∈A

(15)

is the distance of one point from a set. Hence, we calculate the distance between two limit sets, namely the limit set obtained from CSTR simulation and its approximation by POD. Table 2 reports the values of the Hausdorff distance computed for 22, 23 and 24 modes, whereas Fig. 11 reports corresponding phase portraits of the limit cycles referred to the state variables at the outlet. As it is seen from Table 2, the distance monotonically Table 2 Hausdorff distance between the limit set of the CSTR solution and the limit sets of three POD model solutions

Fig. 9. Error at the outlet for conversion degree () and temperature () as a function of the number of modes, Da = 0.035.

No. of modes

Hausdorff distance

22 23 24

9.0324E−003 6.9190E−003 5.4404E−003

K. Bizon et al. / Computers and Chemical Engineering 32 (2008) 1305–1315

1313

Fig. 12. Solution diagram for conversion degree vs. Damk¨ohler. Pe = 100. Stable stationary solutions and stable oscillations (—); unstable oscillations (– – –); unstable stationary solutions (· · ·). Diamonds: POD/Galerkin, 23 modes, snapshots taken for dynamic regimes only.

Fig. 11. Comparison of phase portraits from the 22-mode (a), 23-mode (b) and 24-mode models (c); POD/Galerkin (—); CSTR model (– – –).

decreases as the number of modes is increased. This is reflected by the fact that the limit cycles approach the ideal solution better and better, especially in the relaxation phase (the rightmost segment of the limit cycle). It is here noted that this measure of quality is “static” in that it does not take into account time dependence of the solutions and thus it conceals phase and frequency errors discussed above. However, the evaluation criterion based on distance between limit sets is the most appropriate if we aim at reproducing asymptotic behaviour in terms of phase portraits and bifurcations.

30, however not very accurate. Thus, in order to capture both the catastrophic saddle-node bifurcation at Da = 0.0216 and the Hopf bifurcation at Da = 0.0946 with a not too large number of modes, the basis function set was recalculated by incorporating some snapshots obtained for an extended range of the Damk¨ohler number, namely from steady states. Fig. 13 reports the solution diagram computed for the system by brute force with 23 POD basis functions obtained with this extended sampling. It is seen that the supercritical Hopf bifurcation is captured, although slightly shifted along the bifurcation parameter axis. The method was also tested for two different values of Pe, namely Pe = 50 and 180. For the CSTR cascade approximation, the number of tanks in the series was chosen 25 and 90, corresponding to 50 and 180 ODE, respectively. Fig. 14 reports the comparison of POD/Galerkin with the CSTR cascade model for Pe = 50. It is seen that the dynamical behaviour of the system

4.2. Use of POD/Galerkin for bifurcation analysis We attempt now to reproduce the solution diagram of Fig. 1a. Fig. 12 reports the solution diagram computed for the system by brute force with POD/Galerkin, compared to the diagram obtained by automatic parameter continuation on the CSTR cascade model and already reported in Fig. 1a. It is seen that the results match in almost the whole range of stable limit cycles. The first bifurcation is also captured, whereas the supercritical Hopf bifurcation is not reproduced. It must be recalled that the POD basis functions were obtained by collecting snapshots from stable oscillations only, in the range of Da in between the values 0.03–0.08. A better reproduction of the solution diagram (not reported in here), including a supercritical Hopf around Da = 0.095, was obtained for a number of modes greater than

Fig. 13. Solution diagram for conversion degree vs. Damk¨ohler. Pe = 100. Stable stationary solutions and stable oscillations (—); unstable oscillations (– – –); unstable stationary solutions (· · ·). Diamonds: POD/Galerkin, 23 modes, snapshots taken for dynamic and steady-state regimes.

1314

K. Bizon et al. / Computers and Chemical Engineering 32 (2008) 1305–1315

5. Conclusions

Fig. 14. Solution diagram for conversion degree vs. Damk¨ohler. Pe = 50. Stable stationary solutions and stable oscillations (—); unstable oscillations (– – –); unstable stationary solutions (· · ·). Diamonds: POD/Galerkin, 11 modes.

is simpler than in the case of Pe = 100, and that POD/Galerkin provides a reasonable approximation with as little as 11 modes. On the contrary, 42 modes do not provide acceptable results for Pe = 180 (Fig. 15). As a first comment, we notice that the dynamical scenario is much more complex that that observed at lower values of the Pe number, with presence of multiple stable limit cycles and appearance of NS bifurcations. POD/Galerkin fails to accurately locate the bifurcation value for the supercritical Hopf and cannot capture the higher amplitude periodic branches to their actual extent. The decrease in performance of the reduced model as the Pe number increases can be explained by observing that, for higher values of Pe, the dynamical system progressively loses its dissipative character, and thus the dimension of the global attractor may increase considerably (Robinson, 2001). This results in the need for higher order reduced models.

Fig. 15. Solution diagram for conversion degree vs. Damk¨ohler. Pe = 180. Stable stationary solutions and stable oscillations (—); unstable oscillations (– – –); unstable stationary solutions (· · ·). Diamonds: POD/Galerkin, 42 modes.

In the attempt of building an accurate reduced dynamical model of a pseudohomogeneous tubular reactor with mass recycle, a POD/Galerkin approach was developed and applied to simulations in the oscillatory regimes. Three different ways of comparing the solutions were employed to highlight features of the reduced model. It was found that the cumulative correlation energy is not a reliable criterion of selection of the truncation order for such a system. Due to frequency drift that may arise in the reduced-model solution with respect to the corresponding reference model solution, a conventional norm computed for the difference between time series is non monotonic, thus inadequate to define the truncation order. The Hausdorff distance between limit sets is found to be the best means to choose the truncation order of the reduced model, although this measure does not carry information about time series. POD/Galerkin performance was then tested for three different values of the P´eclet number. For the chosen system, the method does not provide striking results: in general, one-third to one-half modes, with respect to the order of the original system, are needed to provide a reasonable approximation, and computations-per-variable are much heavier that those required for a CSTR cascade model. It was also observed that the method becomes less effective as the P´eclet number increases. In general, the performance increases by incorporating samples from a wider interval of the bifurcation parameter Da and by including all types of solutions.

References Adrover, A., Cerbelli, S., & Giona, M. (2002). A spectral approach to reaction/diffusion kinetics in chaotic flows. Computers and Chemical Engineering, 26(1), 125–139. Adrover, A., & Giona, M. (2003). Modal reduction of PDE models by means of Snapshot Archetypes. Physica D, 182(1), 23–45. Alonso, A. A., Kevrekidis, I. G., Banga, J. R., & Frouzakis, C. F. (2004). Optimal sensor location and reduced order observer design for distributed process systems. Computer and Chemical Engineering, 28, 27–35. Aris, R., Aronson, D. G., & Swinney, H. L. (Eds.). (1991). Patterns and dynamics in reactive media. New York, USA: Springer-Verlag. Aubry, N., Lian, W.-Y., & Titi, E. S. (1993). Preserving symmetries in the proper orthogonal decomposition. SIAM Journal of Scientific Computation, 14, 483–505. Berezowski, M. (2000). Spatio-temporal chaos in tubular chemical reactors with the recycle of mass. Chaos Solitons and Fractals, 11, 1197–1204. Bretsznajder, S., Kawecki, W., & Marcinkowski, S. (1973). Podstawy og´olne technologii chemicznej. Warszawa: Wydawnictwa Naukowo-Techniczne. (in Polish). Bui-Thanh, T., Damodaran, M., & Willcox, K. (2003). Proper orthogonal decomposition extensions for parametric applications in transonic aerodynamics. AIAA Paper 2003-4213. Cizmas, P. G., Palacios, A., O’Brien, T., & Syamlal, M. (2003). Proper orthogonal decomposition of spatio-temporal patterns in fluidized beds. Chemical Engineering Science, 58, 4417–4427. Crommelin, D. T., & Majda, A. J. (2004). Strategies for model reduction: Comparing different optimal bases. Journal of the Atmospheric Sciences, 61, 2206–2217. Danby, S. J., & Echekki, T. (2005). Proper orthogonal decomposition analysis of autoignition simulation data of nonhomogeneous hydrogen–air mixtures. Combustion and Flame, 144, 126–138.

K. Bizon et al. / Computers and Chemical Engineering 32 (2008) 1305–1315 Doedel E. J., Paffenroth R. C., Champneys A. R., Fairgrieve T. F., Kuznetsov Y. A., Sandstede B., et al. (2001). Auto 2000: Continuation and bifurcation software for ordinary differential equations (with HomeCont). Technical report. Druault, P., Guibert, P., & Alizon, F. (2005). Use of proper orthogonal decomposition for time interpolation from PIV data. Experiments in Fluids, 39, 1009–1023. Eriksson, P., Jim´enez, C., B¨uhler, S., & Murtagh, D. (2002). A Hotelling transformation approach for rapid inversion of atmospheric spectra. Journal of Quantitative Spectroscopy and Radiative Transfer, 73, 529–543. Gokulakrishnan, P., Lawrence, A. D., McLellan, P. J., & Grandmaison, E. W. (2006). A functional-PCA approach for analyzing and reducing complex chemical mechanism. Computers and Chemical Engineering, 30, 1093–1101. Graham, M. D., & Kevrekidis, I. G. (1996). Alternative approaches to the Karhunen–Lo`eve decomposition for model reduction and data analysis. Computers and Chemical Engineering, 20, 495–506. Grossman, R. L., Kamath, C., & Kegelmeyer, P. (2001). Data mining for scientific and engineering applications. Dordrecht, NL: Kluwer Academic Publishers. Holmes, P., & Lumley, J. L. (1996). Turbulence, coherent structures, dynamical systems and symmetry (p. 86). Cambridge, UK: Cambridge University Press. Horn, R. A. (2002). Normal matrices with a dominant eigenvalue and an eigenvector with no zero entries. Linear Algebra and its Applications, 357, 35–44. Jolliffe, I. T. (2002). Principal Component Analysis. New York: Springer. Joyner, M. L. (2004). Comparison of two techniques for implementing the proper orthogonal decomposition method in damage detection problems. Mathematical and Computer Modelling, 40, 553–571. Kirby, M., & Sirovich, L. (1990). Application of the Karhunen–Lo`eve procedure for the characterization of human faces. IEEE Transactions on Pattern Analysis Machine Intelligence, 12, 103–108. Kuznetsov, Y. A. (1998). Elements of applied bifurcation theory (2nd ed.). New York, NY: Springer. Levenspiel, O. (1998). Chemical reaction engineering (3rd ed., p. 321). Wiley. Liang, Y. C. (2002). Proper orthogonal decomposition and its applications—Part 2: Model reduction for MEMS dynamical analysis. Journal of Sound and Vibrations, 256, 515–532. Liu, X., & Jacobsen, E. W. (2004). On the use of reduced order models in bifurcation analysis of distributed parameter systems. Computers and Chemical Engineering, 28, 161–169. Lumley, J. L. (1967). The structure of inhomogeneous turbulent flow. In A. M. Yaglom & V. I. Tatarski (Eds.), Atmospheric turbulence and radio wave propagation (p. 160). Moscow: Nauko. Oran, E. S., & Boris, J. P. (1987). Numerical simulation of reactive flow. New York, USA: Elsevier.

1315

Padhi, R., & Balakrishnan, S. N. (2003). Proper orthogonal decomposition based optimal neurocontrol synthesis of a chemical reactor process using approximate dynamic programming. Neutral Networks, 16, 719–728. Phillipson, P. E., & Schuster, P. (2001). Dynamics of relaxation oscillations. International Journal of Bifurcation and Chaos, 11, 1471–1482. Quarteroni, A., & Valli, A. (1994). Numerical approximation of partial differential equations. Springer series in computational mathematics Berlin: Springer. Robinson, J. C. (2001). Infinite-dimensional dynamical systems. Cambridge University Press. Rowley, C. W., Colonius, T., & Murray, R. M. (2001). Dynamical models for control of cavity oscillations. AIAA Paper 2001-2126. Russo, L., Mancusi, E., Continillo, G., & Crescitelli, S. (2001). Reduced basins of attraction of large dynamical systems for the analysis of chemical reactor models. In 5th Italian conference on chemical and process engineering. Shvartsman, S. Y., Theodoropoulos, C., Rico-Martinez, R., Kevrekidis, I. G., Titi, E. S., & Mountziaris, T. J. (2000). Order reduction for nonlinear dynamic models of distributed reacting systems. Journal of Process Control, 10, 177–184. Sirsup, S., Karniadakis, G. E., Xiu, D., & Kevrekidis, I. G. (2005). Equationfree/Galerkin-free POD-assisted computation of incompressible flows. Journal of Computational Physics, 207, 568–587. Smuła, J. (2006). Analysis of relaxation oscillations of tubular chemical reactors with external feedback by selected decomposition methods. Ph.D. thesis. Poland: Silesian University of Technology at Gliwice. Stuart, A. M., & Humphries, A. R. (1998). Dynamical systems and numerical analysis (p. 645). Cambridge University Press. Subramanian, S., & Balakotaiah, V. (1997). Classification of steady-state and dynamic behaviour of a well-mixed heterogeneous reactor model. Chemical Engineering Science, 52, 961–978. Tarazaga, P., Raydan, M., & Hurman, A. (2001). Perron–Frobenius theorem for matrices with some negative entries. Linear Algebra and its Applications, 328, 57–68. Utturkar, Y., Zhang, B., & Shyy, W. (2004). Reduced-order description of fluid flow with moving boundaries by proper orthogonal decomposition. Heat and Fluid Flow, 26, 276–288. Villadsen, J., & Michelsen, M. L. (1978). Solution of differential equation models by polynomial approximation. Englewood Cliffs, NJ: Prentice Hall. Yuan, T., Cizmas, P. G., & O’Brien, T. (2005). A reduced-order model for a bubbling fluidized bed based on proper orthogonal decomposition. Computers and Chemical Engineering, 30, 243–259. Zhang, Y., Henson, M. A., & Kevrekidis, I. G. (2004). Nonlinear model reduction for dynamic analysis of cell population models. Chemical Engineering Science, 58, 429–445.