Ind. Eng. Chem. Res. 2008, 47, 4797–4806
4797
Dynamics Of Axially Dispersed Reactors Mingheng Li* Department of Chemical and Materials Engineering, California State Polytechnic UniVersity, Pomona, California 91768
Axially dispersed reactors with simultaneous convection, dispersion, and reaction might be used to approximate real reactors. While the steady-state solution to an axially dispersed reactor is well-described in the literature, its dynamic behavior is not often discussed. Studying this problem could provide insights into a class of industrial concentration transition problems, such as the product change in a glass melting furnace where one glass product is gradually replaced by another with a different composition. In this work, the dynamics in a tubular reactor with significant axial dispersion are explored following a state-space approach. The statespace models can be constructed either analytically or numerically. In the analytical approach, a continuoustime state-space model is derived from eigenfunction expansion followed by solving a Sturm-Liouville problem. In the numerical approach, a high-dimensional discrete-time state-space model is formulated by the discretization of the process PDE using the implicit finite-difference scheme. The proper orthogonal decomposition (POD) technique is then used to reduce the order of the state-space model. In either case, an optimal trajectory of the feed composition with bounded constraints to achieve a transition in the outlet concentration is solved based on the reduced models using quadratic programming. Computer simulations are used to demonstrate that substantial dispersions, if they exist in a reactor, can significantly affect the reactor dynamics (i.e., a long time might be required to make a concentration transition). In such a case, an optimal feed trajectory solved by dynamic optimization can be used to significantly reduce the transition time. Potential applications to the glass product transition are discussed. 1. Introduction Because it is common for modern chemical plants to produce a series of products with different compositions using the same reactor, the product transition becomes a very important subject at the interface of reactor engineering and optimal control. Representative industrial examples include the transition of glass products of different colors in a glass melting furnace1,2 and the transition of polymer products of different grades in a polyethylene plant.3–5 In a glass melting furnace, the transition time can be up to several days or even weeks, primarily due to the long residence time and the high dispersion rate of the glass melt.1,6 A reduction in the transition time implies more saleable products and less energy consumption but generally relies on an in-depth understanding of the reactor dynamics.7–9 The readers may refer to open literature for recent studies in the area of analysis, control, and optimization of reactor dynamics10–17 and for solution techniques to dynamic optimization of distributed parameter systems using collocation and shooting formulation methods, etc.18–21 This work focuses on the dynamics of an axially dispersed reactor, which can sometimes be used to approximate real reactors.22,23 While the steady-state solution to an axially dispersed reactor is well-described in the literature,24 its dynamic behavior is not often discussed (noticeable exceptions are the dynamic analysis of distributed parameter tubular reactors25 and the transient analysis of axially dispersed reactors with various boundary conditions26). Studying such a problem (generally described by partial differential equations (PDEs) with timevarying boundary conditions) might provide insights into a class of industrial concentration transition problems in which the grade of the final product is determined by the concentration of a key component, which, in turn, is controlled by the feed composition * Tel.: +1-909-869-3668.
[email protected].
Fax:
+1-909-869-6920.
E-mail:
at the entrance of the reactor. The product transition in a glass melting furnace falls into this category. In this case, the key component is typically a colorant. For example, Solargreen and Solextra are two automotive glass products offered by PPG Industries, Inc.1 Because of a difference in the percentage of iron oxide, these products have different aesthetic appearances. Besides iron oxide, other colorant agents such as nickel oxide, cobalt oxide, chromium, and selenium may also be introduced to the glass melt as body tints to adjust the solar performance of a glass product.27 Because of the infinite dimensional nature of the PDE systems that describe the dynamics of such processes, model reduction techniques are generally necessary to convert the original PDE systems to finite dimensional ones described by ordinary differential equations (ODEs) or algebraic equations with reasonable accuracy.18,28–32 In this work, the state-space approach is used and the dynamic models are constructed both analytically and numerically. In the analytical approach, a continuous-time state-space model is derived from eigenfunction expansion followed by solving a Sturm-Liouville problem. The eigenmodes of the PDE system are then divided into a finite dimensional slow set and an infinite dimensional fast complement. Finally, a finite dimensional set of ODEs written in the state-space form is derived to capture the dominant dynamics of the PDE system. In the numerical approach, a highdimensional discrete-time state-space model is formulated by the discretization of the process PDE using the implicit finite difference scheme. The proper orthogonal decomposition (POD) technique is then used to extract the most important characteristic spatial profiles in order to reduce the order of the statespace model. In either case, an optimal trajectory of the feed concentration with bounded constraints to achieve a transition in the outlet concentration is solved based on the reduced models using quadratic programming. Linear quadratic regulator (LQR) theory may also be used in the absence of input
10.1021/ie800083e CCC: $40.75 2008 American Chemical Society Published on Web 06/24/2008
4798 Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008
constraints. Computer simulations are used to demonstrate that substantial dispersions, if they exist in a reactor, can significantly affect the reactor dynamics (i.e., a long time might be required to make a concentration transition). In such a case, an optimal feed trajectory solved by dynamic optimization can be used to significantly reduce the transition time. Implications to the glass product change process are discussed at the end of this work. 2. Problem Description. Consider an isothermal tubular reactor with substantial axial dispersion and a first-order reaction A ⇒ B; the evolution of the concentration of component A is described by a PDE subject to the so-called Danckwerts boundary conditions:24 ∂U(z, t) ) ∂t
2
∂ U(z, t) ∂U(z, t) +D - kU(z, t) ∂z ∂z2 ∂U(z, t) (1) VU(0-, t) ) VU(0+, t) - D s.t. ∂z z)0+ ∂U(z, t) )0 ∂z z)L where t is the time, V is the fluid velocity in the reactor, L is the length of the reactor, D is the dispersion coefficient, U(0-, t) is the inlet concentration, and U(L, t) is the outlet concentration. In this case, the reactor becomes an ideal plug-flow reactor (PFR) and an ideal continuous stirred tank reactor (CSTR) when D approaches zero and infinity, respectively.22 Such a PDE may be written in a dimensionless form as follows: -V
|
|
∂U(z, t) 1 ∂2U(z, t) ∂U(z, t) + - DaU(z, t) ) Pe ∂z2 ∂z ∂t 1 ∂U(z, t) (2) U(0-, t) ) U(0+, t) s.t. Pe ∂z z¯ )0+ ∂U(z, t) )0 ∂z z¯ )1 where τ ) L/V (the characteristic time of the reactor), Pe ) LV/D (the Peclet number), Da ) kL/V (the Damkohler number), jz ) z/L,tj ) t/τ. In the remainder, we focus on the dimensionless form of this dynamic process and drop the bar sign on the variables. The optimal control problem is formulated in conjunction with the glass product transition, and it is applicable to a class of problems that fall into the same category. In this case, the input variable is the colorant concentration in the batch material (or u(t) ) U(0-, t)) and the output variable is the colorant concentration in the final glass product (or y(t) ) U(1, t)). The objective is to penalize a quadratic control energy functional subject to the process dynamics with bounded input constraints. Note that, because of the infinite dimensional nature of such a dynamic process, an exact measure of the transition time is not possible. Specifically, the optimal control problem takes the following form:
|
|
min J ) u(t)
s.t.
∫
∞
0
∂UΩ(z, t) ) ∂t
2 ∂UΩ(z, t) 1 ∂ UΩ(z, t) + ∂z Pe ∂z2
du(t) dt ∂UΩ(z, t) 1 ∂UΩ(z, t) )0 UΩ(0, t) ) Pe ∂z ∂z z)1 z)0+ DaUΩ(z, t) - Dau(t) -
|
s.t.
|
(4) The corresponding eigenvalue problem is 1 φ ′′ (z) - Daφ(z) Pe 1 (5) φ(0) ) φ ′ (0) s.t. Pe φ ′ (1) ) 0 The solution to an eigenvalue problem slightly different from eq 5 (i.e., without the reaction term) has been provided in a book.33 In order to solve eq 5, one can first multiply both sides of the differential equation with a factor σ(z) ) exp(-Pe · z), or -λφ(z) ) -φ ′ (z) +
′
-λ exp(-Pe · z)φ(z) )
[ Pe1 exp(-Pe · z)φ ′ (z)] Da exp(-Pe · z)φ(z) 1 φ(0) ) φ ′ (0) Pe φ ′ (1) ) 0
s.t.
(6) Then the eigenfunctions (not normalized) satisfying the above Sturm-Liouville problem can be derived as φi(z) ) exp 1 Pe · z 2 γi cos(θiz) + sin(θiz) 2 Pe
(
)[
]
(7)
θi > 0
(8)
where θi satisfies (4θi2 - Pe2) sin(θi) ) 4Peθi cos(θi), and the eigenvalues are related to θi as λi ) Da + (θi2 ⁄ Pe) + (Pe ⁄ 4)
(y˜(t))2 + ε2(u˜(t))2dt
∂U(z, t) 1 ∂2U(z, t) ∂U(z, t) )+ - DaU(z, t) ∂t ∂z Pe ∂z2 ∂U(z, t) 1 ∂U(z, t) )0 U(0-, t) ) U(0+, t) Pe ∂z z)0+ ∂z z)1
|
at the entrance might cause imperfect mixing, which leads to quality problems. Moreover, when the optimization is solved as a purely mathematical problem, negative concentrations might be in the optimal solution although they are not possible in practice. Bounded constraints in the inlet concentration will avoid these issues. 3. State-Space Approach Based on Analytical Eigenfunctions. On the basis of the concept of eigenfunction expansion, if U(z, t) ) UΩ(z, t) + u(t), where UΩ(z, t) ) Σi∞) 1 ai(t)φi(z), it can be derived from the original PDE that
|
(9)
According to the Sturm-Liouville theorem, the eigenfunctions satisfying eq 5 are orthogonal to each other with a weight function σ(z) ) exp(-Pe · z),34 and it is derived that (3)
where y˜(t) ) y(t) -yf, u˜(t) ) u(t) - uf (uf and yf are the steadystate concentrations of the key component at the inlet and outlet of the reactor), and ε represents the weight on the control action. Bounded input constraints are generally necessary in the formulation of the optimization problem. For example, in the glass transition process, too high a concentration of the colorant
2 2 2 1 + + θ 〈φi(z), φi(z)〉σ(z) ) 2 Pe Pe2 i 〈φi(z), φj(z)〉σ(z) ) 0 (i * j) 8θi 〈φi(z), 1〉σ(z) ) 2 Pe + 4θi2
(10)
Taking the inner product of φi/(z) with both sides of eq 4 and dividing by 〈φi(z), φi(z)〉σ(z) yields
Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008 4799
i ) 1, 2, · · · , ∞ (11)
˙ ai(t) ) -λiai(t) + bi[-Da · u(t) - ˙ u(t)],
or y ) Gu + Hx(t0), where y ) [y(t1) u(t2)... y(tn)]T], u ) [u(t0) u(t1)... u(tn-1)]T,
where bi ) 〈φi(z), 1〉σ(z)/〈φi(z), φi(z)〉σ(z) The output is as follows: ∞
y(t) ) U(1, t) )
∑
φi(1)ai(t) + u(t)
(12)
i)1
Equation 11 has a derivative of the input variable u(t) that cannot be effectively handled by the standard state-space formulation. To resolve this issue, a new variable xi(t) ) ai(t) + biu(t) is introduced to convert eqs 11 and 12 into the following form: ˙xi(t) ) -λixi(t) + bi(λi - Da)u(t), ∞
∑
y(t) )
i)1
∑
[
cTd Adbd l cTd An-1 d bd
biφi(1)]u(t)
x(tk+1) ) Adx(tk) + bdu(tk) y(tk) ) cTd x(tk) + ddu(tk)
(16)
where Ad and bd are functions of A, b, and a given time interval ∆t.35 cd ) c and dd ) d. As a result, the output at each time step is expressed as a function of the input and the initial state as follows:
[ ][
l
···
···
··
·
··
·
l cTd An-1 d bd
·
··
0 l
·
··
cTd bd
0
··
·
l
·
·· ··
··
cTd bd
· · · cTd Adbd
] [] cTd Ad
cTd A2d l , H) l , l dd cTd And cTd bd
and n∆t is the horizon length. Therefore, the discretized form of the optimal control problem described by eq 3
∑ (y˜(t ))
2
i
+ ε2(u˜(ti-1))2
s.t.
(15) ATS + SA - (Sb + n)R-1(Sb + n)T + Q ) 0 where the sign ∼ represents the nominal value of the variable from its steady-state, Q ) ccT, n ) cd, and R ) d2 + ε2 to be consistent with the optimization functional in eq 3. In the presence of input constraints, the optimal control problem can only be solved numerically. An open-loop control strategy is employed in this work. First, we derive a discretized state-space model:
···
···
i)1
i)1
dd T cd bd · ··
·
∞
x˙(t) ) Ax(t) + bu(t) (14) y(t) ) cTx(t) + du(t) On the basis of the state-space model represented by eq 14, the solution to the optimal control problem described in eq 3 in the absence of input constraints is given by the state feedback law: u˜ ) -kTx˜, where k ) R-1(Sb + n), and S is determined by the Riccati equation:35
cTd Adbd
···
min J ) u˜(t)
As will be shown later, the magnitude of λi increases following an approximately quadratic function. Given the infinite dimensional nature of the state-space system described by eq 13 and the fact that its dynamics are dominated by the slow eigenmodes, it is possible to approximate the original system with a finite dimensional one. Let x(t) ) [x1(t) x2(t)... xr(t)]T, A ) diag{-λ1 -λ2... -λr}, b ) [b1(λ1 - Da) b2(λ2 - Da)... br(λr - Da)]T, c ) [φ1(1) φ2(1)... φr(1)]T, and d ) 1 - Σir ) 1 biφi(1); then the following continuous-time state-space model can be constructed:
cTd bd
···
·
l
(13)
y(t1) y(t2) ) l l y(tn)
dd T cd bd · ··
i ) 1, 2, · · · , ∞
∞
φi(1)xi(t) + [1 -
G)
cTd bd
l dd
]
×
[ ][ ]
· · · cTd Adbd cTd bd
cTd Ad u(t0) u(t1) cTd A2d + l l x(t0) l l u(tn-1) cTd And
(17)
(18)
umin e u e umax is equivalent to the following quadratic programming problem in the form of Euclidean norms: min J ) ||Gu ˜ + Hx ˜(t0)||2 + ε2||u ˜||2 u ˜ (19)
s.t. umin e u e umax
provided that a steady state is reached at the tf ) n∆t.36 From eq 19, it can also be derived that the optimal control trajectory is ˜(t0) ˜ ) -(GTG + ε2I)-1GTHx u
(20)
if there are no input constraints. The stability of the system under the proposed open-loop optimal control should be guaranteed. First of all, the openloop system is stable. Because the input variable solved by the optimization problem is bounded and converges to its new steady state, the output should also be bounded and should converge to its new steady state. Also note that the modes of the infinite dimensional system neglected for controller design are “more” stable; the stability of the original infinite PDE system would not be affected when the controller is applied. A better manipulation of the state variable can be handled by including state constraints in the optimization problem. Interested readers are directed to open literature for model predictive control of PDEs with both state and input constraints.11 4. State-Space Approach Based on Empirical Eigenfunctions. The state-space model of a dynamic PDE system can also be derived using numerical discretization methods such as finite volume/finite difference and finite element, etc. The implicit forward-time, center-space finite difference scheme is used in this work to enhance numerical stability. As a result of the discretization in both time and space of the original PDE, the following algebraic equations are obtained: U(zi+1, tj+1) - U(zi-1, tj+1) U(zi, tj+1) - U(zi, tj) )+ ∆t 2∆z 1 U(zi+1, tj+1) - 2U(zi, tj+1) + U(zi-1, tj+1) Pe (∆z)2 Da · U(zi, tj+1), i ) 1, ..., N + 1 1 U(z2, tj) - U(z0, tj) u ) U(z1, tj) Pe 2∆z U(zN+2, tj) - U(zN, tj) 0) 2∆z
(21)
where N + 1 is the total number of nodes in the spatial domain (with nodes 1 and N + 1 representing the left and right
4800 Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008
boundaries), and i and j are the spatial and temporal indices, respectively. Note that two fictitious points (z0 and zN+2) are used to approximate ∂U/∂z on the boundaries. The values of the variable U at these fictitious points can be expressed through boundary conditions as U(z0, tj) ) -γU(z1, tj) + U(z2, tj) + γu U(zN, tj) U(zN+2, tj) )
(22)
where γ ) 2Pe(∆z) (note that the PDE is in the dimensionless form, and therefore, ∆z is dimensionless). Sorting and combining terms in the above equation yields: U(zi, tj) ) (-R - β)U(zi-1, tj+1) + (1 + 2β + κ)U(zi, tj+1) + (R - β)U(zi+1, tj+1), i ) 2, ..., N U(z1, tj) ) [1 + 2β + κ + (R + β)γ]U(z1, tj+1) - 2βU(z2, tj+1) [(R + β)γ]u U(zN+1, tj) ) -2βU(zN, tj+1) + (1 + 2β + κ)U(zN+1, tj) (23) where R ) ∆t/2∆z, β ) ∆t/Pe(∆z)2, and κ ) Da∆t. Let x(k) ) [U(z1, tk) U(z2, tk)... U(zN+1, tk)]T, b ) [(R + β)γ 0... 0]T, c ) [0 0... 0]T,
[
An ) k1 -2β -R - β 1 + 2β + κ · ··
R-β · · ·· ·· -R - β 1 + 2β + κ R-β -2β 1 + 2β + κ
]
where k1 ) 1 + 2β + κ + (R + β)γ; the original PDE can be converted to a state-space form of the following: Anx(k + 1) ) Aox(k) + bu(k) y(k) ) cTx(k)
(24)
The resulting reduced model is as follows,
where Ao ) I, which is the identity matrix (note Ao might be a sparse matrix similar to An if the semi-implicit numerical scheme is used). Because An is high-dimensional, a large amount of computational cost might be involved in the optimization. To circumvent the computational complexity brought by this, the Galerkin projection is applied on the high-dimensional statespace model.37 First, we define the ensemble of U(z, t) using the snapshots generated by the open-loop simulation of the highdimensional state-space model:38
[
U(z1, t1) U(z1, t2) · · · U(z1, tK) U(z2, t2) · · · U(z2, tK) U(z2, t1) U) l l U(zN+1, t1) U(zN+1, t2) · · · U(zN+1, tK)
Figure 1. First 50 eigenvalues and first 5 normalized eigenfunctions derived from the eigenfunction expansion.
]
a(k + 1) ) Ara(k) + bru(k) crTa(k) y(k) )
where Ar ) (ΦTAnΦ)-1ΦTAoΦ, br ) (ΦTAnΦ)-1ΦTb, and crT ) cTΦ, provided that ΦTAnΦ (or just An) is full rank. Therefore, the order of the state-space model is reduced from N + 1 to r through the Galerkin projection. The reduced model can then used for control and optimization purposes. The larger the r is, the more ensemble information of the snapshots is retained in the model, however, the more computational cost is required. In practice, r is chosen so that (|x| - |Φa|)/|x| is less than a user-specified value. The optimal control problem based on the reduced model described by eq 26 is the following, ∞
A singular value decomposition (SVD) of U yields that U ) 1,i)j Σi ) 1 φiλiψiT, where φiTφj ) {0,i+j . Let Φ (Φ ) [φ1, φ2,..., φr]) be the matrix composed by empirical eigenfunctions (φi) corresponding to the first r singular values of the matrix composed by the ensemble of U; it can be shown that, within the generated data set, error ) (|x| - |Φa|)/|x| ) 1 - Σir ) 1 λi/ ΣiK) 1 λi, where the eigenvalues satisfy λ1 g λ2 g... g λr g... g 0. If the snapshots are adequately representative of the process dynamics, one can approximate the variable using the weighted sum of the first reigenfunctions with reasonable accuracy, i.e., x(k) ≈ Φa(k), and a(k) ≈ ΦTx(k). Therefore, eq 24 can be converted to the following form:
min J ) u(t)
K
ΦTAnΦa(k + 1) ) ΦTAoΦa(k) + ΦTbu(k) y(k)
)
T
c Φa(k)
(26)
(25)
∑ (a˜ Qa˜ + u˜ Ru˜) T
T
i)1
s.t.
(27)
umin e u e umax where Q ) crcrT ) ΦTccTΦ and R ) ε2 in order to be consistent with the original optimal problem described in eq 3. Again, the tilde sign represents the deviation of the variable from its steadystate value. If there are no input constraints, the solution to this optimization problem is given by the state feedback law: u˜ ) -kTa˜, where k ) (brTSbr + R)-1(ArTSbr) and S is determined by the discrete Riccati equation:39 ArTSAr - S - (ArTSbr)(brTSbr + R)-1(brTSAr) + Q ) 0 (28)
Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008 4801
In the presence of input constraints, quadratic programming can be used following a similar approach presented in Section 3. Note that the state-space model constructed using finitedifference is already in the discrete-time form and no discretization is necessary. Moreover, dd in this case is zero. Therefore, if one set u(ti) ≡ 1, the step response can be easily found to be
[ ][ y(t1) y(t2) ) l l y(tn)
cTd bd
0
...
cTd Adbd
cTd bd · ··
·
l
·
··
l ···
0
··
·
cTd An-1 d bd
...
··
···
l ·
·· T cd bd T cd Adbd
l 0 cTd bd
][ ] [ ] P(t1) 1 P(t2) 1 l ) l l l P(tn) 1
(29)
where P(ti) is the step response of the system at t ) ti; it can then be determined that
[ ][ ] cTd bd T cd Adbd l l
or
G)
[
cTd An-1 d bd
P(t1) P(t2) - P(t1) l l P(tn) - P(tn-1)
P(t1) P(t2) - P(t1) ) l l P(tn) - P(tn-1)
0 ··· P(t1) ··· ·
··
·
··
·
···
(30)
···
0 l
·
l
··
P(t1) ·· 0 · · · P(t2) - P(t1) P(t1)
]
(31)
The matrix G has exactly the same form as the matrix ∆P developed in an input/output approach,36 which suggests that these two approaches are inherently consistent with each other. However, the information of the state is only available in the state-space model, and therefore, it is more flexible than the input/output model.36 5. Results and Discussion The parameters used in the simulations are Pe ) 10 and Da ) 1. The eigenvalues and the corresponding eigenfunctions are solved based on the method of eigenfunction expansion presented in Section 3, and the results are shown in Figure 1. Clearly, θi increases with i almost linearly. As a result, λi increases with i following an approximately quadratic function (recall that λi ) Da + (θi2/Pe) + (Pe/4)). Because the system dynamics are mainly dependent on the slow eigenmodes (associated with λi close to zero) as the fast ones flatten out quickly, only a finite number of eigenmodes are required to approximately describe the process dynamics. For the analytical eigenfunctions, it is interesting to observe that φi(z) crosses the horizontal axis (i - 1) times and the maximum of φi(z) occurs at z ) 1. Moreover, the larger the i, the larger is the φi(1). A rigorous mathematical analysis might help explain this, but it is out of the scope of this work. When the implicit finite difference is used, 2001 snapshots are generated from a step response of the dynamic system. The SVD is then applied to the ensemble of the state variable x. The singular values are normalized and sorted in a
Figure 2. First 20 eigenvalues and first 5 eigenfunctions derived from an ensemble of the concentration profile based on the state-space model derived from implicit finite difference.
descending order. The normalized eigenvalues derived from POD and their cumulated sums (Σmi ) 1 λm/ΣmK ) 1 λm) are shown in Figure 2a. The first five empirical eigenfunctions corresponding to the five largest eigenvalues λi are shown in Figure 2b. It is seen that the process can be described using about 10 empirical eigenfunctions with very reasonable accuracy. However, it is important to note that a reduced model derived from POD can only be as good as the snapshot ensemble information. In order for the reduced model to be accurate for control and optimization purposes, the system has to be adequately excited so that the snapshot ensemble possesses adequate representative information of the process dynamics. To verify this, an analysis of the quality of the derived reduced model is conducted using the Bode plot,40,41 and the results are shown in Figure 3. It is seen that the reduced model derived from POD is consistent with the highdimensional model at low frequencies, while observable differences occur at frequencies higher than 20 Hz. This is acceptable for application purposes because varying the input variable with a very high frequency is unpractical in many industrial processes (e.g., the batch composition in the glass product transition). The reduced model derived from POD is also tested using two input functions, and the responses are compared with those based on the high-dimensional model (order ) 201). As shown in Figure 4, the two models are consistent with each other. It is also noticed that, while all the analytical eigenfunctions used in the eigenfunction expansion share the same boundary conditions at z ) 0 and z ) 1 (see Figure 1), the same behavior is only observed in the empirical eigenfunc-
4802 Ind. Eng. Chem. Res., Vol. 47, No. 14, 2008
Figure 3. Bode plots of the high-dimensional and reduced models derived from eigenfunction expansion and finite difference with POD.
tions derived from finite difference and POD at z ) 1. This is because the variable of interest in the series expansion is U(z, t) - u(t) in the former case but just U(z, t) in the latter, and the boundary condition of U(z, t) at z ) 0 is changing with time, i.e., u(t) ) U(0+, t) - (1/Pe) (∂U(z, t)/∂z)|z ) 0+ (or its equivalent discretized form described by eq 23 based on the finite difference scheme). To understand better the empirical eigenfunctions derived from POD and finite difference, they are just a set of characteristic profiles of the linear system described by eq 24, which in turn, is an “approximation” of the original PDE described by eq 2 with a resolution of ∆t in time and ∆z in space due to the finite discretization nature. Nevertheless, these two formulations based on analytical and empirical eigenfunctions are very close to each other for application purposes in which the frequency is
