Fast Estimation of Quasi-Steady States of Cyclic Nonlinear Processes

266

Ind. Eng. Chem. Res. 2006, 45, 266-291

Fast Estimation of Quasi-Steady States of Cyclic Nonlinear Processes Based on Higher-Order Frequency Response Functions. Case Study: Cyclic Operation of an Adsorption Column Menka Petkovska* and Ana Markovic´ † Department of Chemical Engineering, Faculty of Technology and Metallurgy, UniVersity of Belgrade, KarnegijeVa 4, 11000 Beograd, Serbia and Montenegro

A new method for fast approximate calculation of quasi-steady states of cyclic processes is presented. The method is based on the concept of higher-order frequency response functions. The system input is represented in the form of Fourier series, whereas the output is presented in the form of Volterra series. For practical applications, both the input and the output series are approximated by finite-length sums. In this way, the approximate periodic quasi-steady state of the system output is calculated directly, without long numerical integrations. Cyclic operation of an adsorption column with periodic fluctuations of the inlet concentration or/and adsorbent temperature is used as a case study for testing the new method. The necessary frequency response functions (FRFs), up to the third order, are derived, based on the equilibrium dispersion model. The method is tested for sinusoidal and rectangular input changes. The approximate solutions based on the FRFs, up to the third order, and a finite number of input harmonics, are calculated for different input frequencies and amplitudes and compared with the numerical solutions. Very good agreement is obtained. 1. Introduction 1.1. Cyclic ProcessessSignificance of Fast Estimation of Quasi-Steady States. The majority of processes in chemical and other process industries are operated in a continuous way, in optimal steady states, which is attained using different process control strategies. Nevertheless, some processes, such as adsorption or ion exchange, are usually performed in fixed beds, for which continuous steady-state operations are not possible. These processes are usually performed in a cyclic way, with adsorption and desorption (regeneration) half-cycles. Several different cyclic adsorption processes have been developed in the last 30 years, such as temperature swing adsorption (TSA), pressure swing adsorption (PSA),1 parametric pumping,2 and cyclic zone adsorption (CZA),3 as well as the simulated moving bed (SMB) technique.4 All these processes were designed to simulate continuous operations; however, in their essence, they are cyclic processes. On the other hand, with the development of process dynamics and modeling, an idea about deliberate unsteady operation of processes that could otherwise be performed in steady-state modes became popular among researchers in the 1960s and 1970s. An overview of the research in this field from that period, including cyclic operation of absorption columns, heat exchangers, and chemical reactors, can be found in a classic book by Douglas.5 The essence is that, for some cases, because of the system nolinearity, a periodic (cyclic) operation, caused by forced periodic modulation of one or more input variables, can be superior to the optimal steady-state operation. The research in this area was mainly focused on periodic operations of catalytic chemical reactors. A comprehensive review of a large number of investigations of catalytic reactors with composition modulation was given by Silveston.6 Catalytic reactors with pressure and temperature modulations have also been investi* To whom all correspondence should be addressed. Tel: (+38111) 3303-610. Fax: (+381-11) 3370-387. E-mail: [email protected]. † Current address: Max-Planck Institute for Dynamics of Complex Technical Systems, Sandtostrasse 1, Magdeburg, Germany.

gated,7,8 as well as periodic operations of trickle-bed reactors.9-11 Investigation of the periodic operation of other processes (e.g., distillation12,13) has also been reported. With the increasing demands for higher process productivity and efficiency, combined with the development of modern and efficient control strategies, we could expect faster development and wider practical application of cyclic processes in the near future. To reduce the experimental work in such investigations, mathematical modeling and simulation is often used as an aid in the analysis and optimization of cyclic processes. These models are usually very complex, in the form of coupled nonlinear partial differential equations, and can generally be solved only numerically. To simulate the periodic quasi-steady state (cyclic) behavior of the process, which is of interest, it is necessary to run the numerical simulation until the transient response completely vanishes. For slow processes, these simulations can be very long. Direct calculation of the periodic quasisteady state, without going through the transient solution of the model equations would be very desirable. Some attempts in this direction have been reported by LeVan and co-workers.14-17 In this paper, we propose a new method for the fast approximate calculation of the quasi-steady states of periodic processes, which avoids long numerical solutions of partial differential model equations. The method is based on Volterra series expansions and the concept of higher-order frequency response functions, which will be briefly explained in the next section. 1.2. Nonlinear Frequency Response and the Concept of Higher-Order Frequency Response Functions. Frequency response (FR) is a quasi-stationary response of a stable system to a periodic (sinusoidal or cosinusoidal) input change around a steady state. Nonlinear frequency response is the FR of a nonlinear system. One of the most convenient tools for treating nonlinear FRs is the concept of higher-order frequency response functions (HFRFs),18 which is based on Volterra series and generalized Fourier transform. This concept will be briefly presented below.

10.1021/ie0505965 CCC: $33.50 © 2006 American Chemical Society Published on Web 11/29/2005

Ind. Eng. Chem. Res., Vol. 45, No. 1, 2006 267

Let us consider a stable system with a single input x and a single output y. Dynamic response of a linear system to an arbitrary input x(τ) can be defined using a convolution integral:

y(τ) )

∫-∞∞ g(τ1)x(τ - τ1) dτ1

(1)

where g(τ) is the so-called impulse-response function of the system, or its kernel. By applying Fourier transform to the function g(τ), the frequency-response function (FRF), or frequency transfer function, is obtained:

G(ω) )

∫-∞∞ g(τ)e-jωτ dτ

(2)

This function is directly related to the amplitude and phase of the system’s quasi-stationary response to a single harmonic input:

x ) A cos(ωτ) w τ f ∞: y(τ) ) A|G(ω)| cos(ωτ + arg(G(ω))) (3) On the other hand, the response of a weakly nonlinear system, for which the system nonlinearity has a polynomial form (or can be developed in a Taylor series), can be represented in the form of a Volterra series:18 ∞

y(τ) )

∑yn(τ)

(4)

n)1

with the nth element of the series defined as

yn(τ) )

∫-∞∞ ‚‚‚ ∫-∞∞ gn(τ1,...,τn)x(τ - τ1) ‚‚‚

x(τ - τn) dτ1 ‚‚‚ dτn (5)

where gn(τ1,...,τn) is the nth-order Volterra kernel, or the generalized impulse response function of order n. The first element of the Volterra series, y1, corresponds to the linearized model, while y2, y3, ... are the “correction” functions of the first, second, ... order. Similarly to Taylor or Fourier series expansions, a Volterra series of indefinite length is needed for exact representation of a nonlinear system; however, for practical applications, finite series can be used. By applying multidimensional Fourier transform on the function gn(τ1,...,τn), the nth-order FRF is obtained:

Gn(ω1,...,ωn) )

∫-∞∞ ‚‚‚ ∫-∞∞ gn(τ1,...,τn)ej(ω τ +‚‚‚ω τ ) dτ1 ‚‚‚ dτn 1 1

n n

(6)

This function is directly related to the nth element of the output in its quasi-steady periodic state, presented in the Volterra series form (eq 4). If the input is a periodic function of a more general form,

∑ Ake

In our previous work, we have been using the concept of HFRFs to investigate adsorption kinetics and equilibria.19-24 The most important results of this investigation are that it is possible to discriminate between different kinetic mechanisms reliably using HFRFs22,24 and that both kinetic and equilibrium parameters, including those which define the system nonlinearity, can be estimated from the same set of experimental data.20,21 The investigation presented in this paper is the first attempt to apply this concept to the estimation of quasi-steady states of periodic processes. 2. Calculation of Periodic Quasi-Steady States Based on Higher-Order Frequency Response Functions 2.1. Model Definition for a Weakly Nonlinear System with Multiple Inputs. In performing periodic operations, generally one or more input variables are modulated in a periodic way. To develop a method for fast calculation of cyclic quasi-steady states, we will consider a system with two periodically modulated inputs x and z. A block diagram, representing a nonlinear system with two inputs, is presented in Figure 1. Note that, for this case, to define the complete model, it is necessary to define three sets of HFRFs: two of them relating the output to each of the inputs, and one cross-function, relating the output to both inputs. This third set contains only functions of the second and higher orders. In Figure 1, we denote the FRFs in the following way: Gn,xn and Gn,zn are the nth-order FRFs corresponding to the individual inputs x and z, respectively, whereas Gn,xmzn-m is the nth-order cross-function, with order m regarding the input x and order n-m regarding the input z. This analysis could be extrapolated to nonlinear systems with three or more input variables. For example, for a nonlinear system with three inputs, it would be necessary to define seven sets of HFRFs: three that relate the output to each of the inputs, three sets of cross-functions that relate the output to each combination of two inputs, and one set of cross-functions that relate the output to all three inputs. This last set would contain only the FRFs of the third and higher orders. 2.2. Calculation of the Periodic Quasi-Steady State of the System Output. It is well-known that any periodic function can be represented as an indefinite sum of harmonic functions, i.e., in the Fourier series form. Accordingly, the input functions x and z can be written in the following way: ∞

x(τ) )

K

x(τ) )

Figure 1. Block diagram of a nonlinear system with two inputs and one output.

jωkτ

(7)

∞

k)1

this nth element of the Volterra series becomes K

yn(τ) )

K

∑∑

k1)1 k2)1

K

‚‚‚

∑ Ak Ak

kn)1

1

2

‚‚‚ Akn ×

Gn(ωk1,ωk2,...,ωkn)ej(ωk1+ωk2+‚‚‚+ωkn)τ (8)

z(τ) )

1

Akejw τ ∑ 2 k)-∞ 1

k

Bkeju τ ∑ 2 k)-∞ k

(9)

(10)

where A-k ) conj(Ak), B-k ) conj(Bk), wk ) kw1, and uk ) ku1, whereas w1 and u1 are the basic frequencies of the input changes x and z, respectively. The complex magnitudes of Ak and Bk are the amplitudes and their phase angles represent the

268

Ind. Eng. Chem. Res., Vol. 45, No. 1, 2006

phases of the harmonics of the corresponding frequencies wk and uk of the inputs x and z, respectively. For the case of two periodic inputs, the periodic quasi-steady state of the output can be represented as a sum of three terms, each of them corresponding to the contribution of each set of HFRFs, presented in Figure 1:

y(τ) ) yx(τ) + yz(τ) + yxz(τ)

(11)

Each of these three terms can be represented by an indefinite Volterra series: ∞

yx(τ) )

∑ yx,n(τ)

(12)

n)1 ∞

yz(τ) )

∑ yz,n(τ) n)1

(13)

∞

yxz(τ) )

∑ yxz,n(τ)

(14)

n)1

with the nth elements of these series defined in the following way: ∞

1

∞

∞

Ak Ak ‚‚‚Ak ∑ ∑ ‚‚‚ k ∑ n k )-∞ k )-∞ )-∞

yx,n(τ) )

1

2

1

2

2

n

×

n

Gn,xn(wk1,wk2,...,wkn)ej(wk1+wk2+‚‚‚+wkn)τ (15) yz,n(τ) )

∞

1

∞

∑ ∑

∞

∑

Bk1Bk2‚‚‚Bkn × kn)-∞ 2n k1)-∞ k2)-∞ Gn,zn(uk1,uk2,...,ukn)ej(uk1+uk2+‚‚‚+ukn)τ (16)

yxz,n(τ) )

1

n-1

∞

‚‚‚

∞

Ak ‚‚‚Ak Bk ∑ ∑ ‚‚‚ k ∑ n m)1 k )-∞ )-∞

‚‚‚Bkn × 1 m m+1 2 1 n Gn,xmzn-m(wk1,...,wkm,ukm+1,...,ukn)ej(wk1+‚‚‚+wkm+ukm+1+‚‚‚+ukn)τ (17)

The Volterra series, given by eqs 11-17, with indefinite sums, represent the periodic quasi-stationary state of the output exactly. In practice, only finite sums can be calculated, resulting with equations corresponding to approximate calculation of the periodic output. The inputs are approximated by finite sums of only the first K harmonics: K

xapp(τ) ) zapp(τ) )

∑

1

k)-K

2

K

1

Akejwkτ

(18)

Bkeju τ ∑ 2 k)-K k

(19)

and the output is approximated by the sum of only the first N elements of the Volterra series: N

yapp(τ) )

N

N

∑ yx,n,app(τ) + n)1 ∑ yz,n,app(τ) + n)1 ∑ yxz,n,app(τ) n)1

(20)

with

K

1

yx,n,app(τ) )

2

K

K

Ak Ak ‚‚‚Ak ∑ ∑ ‚‚‚ k ∑ n k )-K k )-K )-K 1

1

2

2

n

×

n

Gn,xn(wk1,wk2,...,wkn)ej(wk1+wk2+‚‚‚+wkn)τ (21) K

1

yz,n,app(τ) )

K

K

Bk Bk ‚‚‚Bk ∑ ∑ ‚‚‚ k ∑ n k )-K k )-K )-K

2

1

1

2

2

n

×

n

Gn,zn(uk1,uk2,...,ukn)ej(uk1+uk2+‚‚‚+ukn)τ (22) yxz,n,app(τ) )

1

n-1

K

∑ ∑

K

‚‚‚

∑

Ak1‚‚‚AkmBkm+1‚‚‚Bkn × kn)-K 2n m)1 k1)-K Gn,xmzn-m(wk1,...,wkm,ukm+1,...,ukn)ej(wk1+‚‚‚+wkm+ukm+1+‚‚‚+ukn)τ (23)

The quality of the approximative solution obtained in this way, i.e., its closeness to the exact solution, increases as both K and N increase. 2.3. Procedure for Approximate Calculation of the Periodic Quasi-Steady State. Based on the analysis presented in the previous two sections, we suggest the following procedure for approximate calculation of the quasi-steady states of periodic processes: Step 1: Postulating of an appropriate mathematical model for the investigated process; Step 2: Derivation of the HFRFs up to the Nth order, based on the postulated model; Step 3: Definition of the periodic input or inputs and their approximation by finite sums, taking into account only the first K harmonics; and Step 4: Calculation of the approximate output, using eqs 2023. In this procedure, the most complex and time-consuming is the second step, the derivation of the HFRFs. Nevertheless, this step is performed only once for each investigated system. Once derived, these functions can be used for any shapes of the periodic input changes and for any set of model parameters. Steps 3 and 4, naturally, must be repeated for each particular simulation. Step 3 can be easily performed by applying the Fourier transform to the periodic input functions, whereas the last step, calculation of the periodic output in its quasi-steady state, practically reduces to simple algebra. To illustrate this procedure and our new method, we will apply it to simulations of the cyclic operation of an adsorption column with periodic modulation of the inlet concentration or/ and the temperature of the entire adsorbent bed. 3. Case Study: Cyclic Operation of an Adsorption Column There were several reasons for choosing an adsorption column as a case study for our new method. First, we were looking for a rather simple system as a first illustration of application of a new method. Second, adsorption columns must be operated periodically, so this application definitely has a practical value. The choice of the inputs that are periodically modulatedsthe inlet concentration and the temperature of the entire adsorbent bedscorrespond to the so-called direct, or standing-wave mode of operation.3 On one hand, the model equations for this case are rather simple, and on the other, it is possible to explore the response of the column to one lumped (the inlet concentration) and one distributed (the temperature of the entire adsorbent bed) input periodic change.


τ e 0, ∀z: C(z) ) Ci(z) ) Cs ) Ci,s, Q(z) ) Qs ) Φ(Cs,Ts), T ) Ts (28) For convenience, the model equations 24, 25, 26, 27, and 28 were transformed to the following nondimensional equations, respectively:

∂c ∂q ∂c 1 ∂2c +f + ) ∂t ∂t ∂x 2NTP ∂x2 x ) 0: c(0,t) ) ci(t) +

∂c | | )0 ∂x |x)1

x ) 1: Figure 2. Three cases of periodic operation of the adsorption column.

The periodic behavior of the adsorption column will be analyzed for the simplest case of a single adsorbing component and for a constant flow rate of the moving phase. Three cases of periodic operation of the adsorption column will be analyzed: Case 1: For periodic modulation of the inlet concentration and constant temperature of the adsorbent bed; Case 2: For periodic modulation of the temperature of the adsorbent bed and constant inlet concentration; and Case 3: For simultaneous periodic modulations of the inlet concentration and the temperature of the adsorbent bed. These three cases are schematically shown in Figure 2. In this figure, Ci is the inlet concentration, Co the outlet concentration, and T the temperature of the bed. Π1(w,t) and Π2(u,t) denote periodic functions of time with frequencies w and u, respectively. The subscript s corresponds to steady-state values. 3.1. Model Equations. As mentioned in Section 2.3. the first step in our procedure for approximate calculation of the periodic quasi-steady state is the formulation of a suitable nonlinear mathematical model. For the adsorption column used as a case study, we will use the equilibrium-dispersion model with a nonlinear adsorption isotherm.25 This model assumes that different mass-transfer resistances and axial dispersion can be represented by a single apparent or effective dispersion coefficient and that local equilibrium between the solid and the fluid phase exists in each point of the column. In our model, we also assume uniform temperature in the entire column. With all these assumptions, the material balance of a single adsorbing component in the adsorption column is obtained in the form of a second-order partial differential equation:

∂C 1 - ∂Q ∂2C ∂C + +V ) Deff 2 ∂τ ∂τ ∂z ∂z

(24)

z ) 0: C(0,τ) ) Ci(τ) + Deff z ) L:

∂C | | )0 ∂z |z)L

∂C | | ∂z |z)0

(25) (26)

and with local equilibrium defined by the adsorption isotherm relation, which is generally nonlinear:

∀z, ∀τ: Q ) Φ(C,T)

(30) (31)

∀x, ∀t: q ) a˜cc + a˜θθ + b˜ ccc2 + b˜ θθθ2 + b˜ cθcθ + c˜cccc3 + c˜θθθθ3 + c˜ccθc2θ + c˜cθθcθ2 + ‚‚‚ (32) t e 0, ∀x: c(x) ) ci(x) ) q(x) ) θ ) 0

(33)

The nondimensional independent variables were defined as

V L

(34a)

z L

(34b)

t)τ x)

the nondimensional concentrations and temperature were defined as

c)

C - Cs Cs

(35a)

q)

Q - Qs Qs

(35b)

ci )

Ci - Cs Cs

(35c)

θ)

T - Ts Ts

(35d)

and the nondimensional parameters (the number of theoretical plates NTP and the porosity factor f) were defined as

NTP ) f)

with the following boundary conditions:

1 ∂c | | 2NTP ∂x |x)0

(29)

VL 2Deff

1 - Qs Cs

(36)

In eq 32, the nondimensional adsorption isotherm relation is presented in the form of a Taylor series expanded around c ) 0 and θ ) 0 (C ) Cs and T ) Ts), with its coefficients defined as

a˜c )

Cs ∂Φ | | ∂C |C)Cs Qs

(37a)

Ts ∂Φ | |C)C ∂T | s Qs

(37b)

T)Ts

(27)

The initial conditions are based on the assumption that the solid and the fluid phases are initially in equilibrium:

a˜θ )

T)Ts

270


b˜ cc )

Cs 1 ∂2Φ | |C)C 2 s 2 ∂C | Qs

(37c)

T)Ts

first-order FRF:

2

Ts 1 ∂2Φ | b˜ θθ ) | 2 ∂T2 |C)Cs Qs

(37d)

G1,c(w1) ) U1(w1)eλ1(w1) + U2(w1)eλ2(w1)

(37e)

where λ1(w1) and λ2(w1) are the characteristic values of the differential equation defining the first-order FRF while U1(w1) and U2(w2) are the integration constants, which can be found in Section A1.1 of the Appendix.

T)Ts

CsTs ∂Φ| b˜ cθ ) |C)C s ∂C∂T | Qs 2

T)Ts

Cs 1 ∂ 3Φ | |C)C 3 s 6 ∂C | Qs

(37f)

T)Ts

G2,cc(w1,w2) ) D1(w1,w2)eµ1(w1,w2) + D2(w1,w2)eµ2(w1,w2) + d1(w1,w2)e(λ1(w1)+λ1(w2)) + d2(w1,w2)e(λ1(w1)+λ2(w2)) +

3

c˜θθθ )

(38)

second-order FRF:

3

c˜ccc )

3.2.1. The Gc Functions. This set of FRFs corresponds to the modulation of the inlet concentration. The following expressions for the first-, second-, and third-order FRFs were obtained:

Ts 1 ∂3Φ | | 6 ∂T3 |C)Cs Qs

(37g)

T)Ts

Cs2Ts 1 ∂3Φ | c˜ccθ ) | 2 ∂C2∂T |C)Cs Qs

(37h)

T)Ts

Here, µ1(w1,w2) and µ2(w1,w2) are the characteristic values corresponding to the differential equation defining the secondorder FRFs and D1(w1,w2), D2(w1,w2), d1(w1,w2), d2(w1,w2), d3(w1,w2), and d4(w1,w2) are the integration constants. Their expressions are all listed in Section A1.2 of the Appendix.

third-order FRF:

2

CsTs 1 ∂3Φ | c˜cθθ ) | C)C 2 ∂C∂T2 | s Qs

d3(w1,w2)e(λ2(w1)+λ1(w2)) + d4(w1,w2)e(λ2(w1)+λ2(w2)) (39)

(37i)

T)Ts

These coefficients can be regarded as local nondimensional first, second, and third derivatives of the adsorption isotherm. 3.2. Adsorption Column Frequency Response Functions. The second step in our procedure for approximate calculation of the periodic quasi-steady state of the adsorption column is derivation of the necessary FRFs. Because we are looking at the periodic behavior of the adsorption column with one or two modulated inputs, we need three series of HFRFs, as explained in Section 2.1. We will denote the nth-order FRFs that correlate the outlet and inlet concentrations with Gn,cn(w1,...,wn), the FRFs that correlate the outlet concentration and the temperature of the adsorbent bed with Gn,θn(u1,...,un), and the mixed FRFs with Gn,cmθn-m(w1,...,wm,um+1,...,un). In the further text, we will call these three sets of FRFs the Gc, Gθ, and Gcθ functions, respectively. The procedure for derivation of HFRFs is rather standard and can be found in our previous papers.21,23 The basic steps of this procedure applied to our case are as follows: (1) Define the inputs, which are the inlet concentration ci(t) and the temperature of the column θ(t), in the form of general periodic functions (such as those in eqs 9 and 10); (2) Represent the concentration in the fluid phase at position x, c(x,t), in the Volterra series form (analogous to eqs 11-17); (3) Express the concentration in the solid phase at position x, q(x,t), using the assumption of local equilibrium and the equilibrium relation (eq 32); (4) Substitute the expressions for ci, θ, c, and q into the model equation (eq 29) and the boundary conditions (eqs 30 and 31); (5) Apply the method of harmonic probing to the equations obtained in step 4 and the solve the resulting equations. The Gc functions (the FRFs with regard to the inlet concentration) have been derived in one of our previous publications.26 Some details about the FRFs derivation can also be found in the Appendix of this manuscript. In the manuscript itself, we will give only the final expressions for the first, second, and third-order FRFs of the adsorption column.

G3,ccc(w1,w2,w3) ) V1(w1,w2,w3)eχ1(w1,w2,w3) + V2(w1,w2,w3)eχ2(w1,w2,w3) + V1(w1,w2,w3)e(µ1(w1,w2)+λ1(w3)) + V2(w1,w2,w3)e(µ2(w1,w2)+λ1(w3)) + V3(w1,w2,w3)e(µ1(w1,w2)+λ2(w3)) + V4(w1,w2,w3)e(µ2(w1,w2)+λ2(w3)) + V5(w1,w2,w3)e(µ1(w1,w3)+λ1(w2)) + V6(w1,w2,w3)e(µ2(w1,w3)+λ1(w2)) + V7(w1,w2,w3)e(µ1(w1,w3)+λ2(w2)) + V8(w1,w2,w3)e(µ2(w1,w3)+λ2(w2)) + V9(w1,w2,w3)e(µ1(w2,w3)+λ1(w1)) + V10(w1,w2,w3)e(µ2(w2,w3)+λ1(w1)) + V11(w1,w2,w3)e(µ1(w2,w3)+λ2(w1)) + V12(w1,w2,w3)e(µ2(w2,w3)+λ2(w1)) + V13(w1,w2,w3)e(λ1(w1)+λ1(w2)+λ1(w3)) + V14(w1,w2,w3)e(λ1(w1)+λ2(w2)+λ1(w3)) + V15(w1,w2,w3)e(λ1(w2)+λ2(w1)+λ1(w3)) + V16(w1,w2,w3)e(λ1(w3)+λ2(w1)+λ2(w2)) + V17(w1,w2,w3)e(λ1(w1)+λ1(w2)+λ2(w3)) + V18(w1,w2,w3)e(λ1(w1)+λ2(w2)+λ2(w3)) + V19(w1,w2,w3)e(λ1(w2)+λ2(w1)+λ2(w3)) + V20(w1,w2,w3)e(λ2(w1)+λ2(w2)+λ2(w3)) (40) where χ1(w1,w2,w3) and χ2(w1,w2,w3) are the characteristic values of the corresponding differential equation and V1(w1,w2,w3), V2(w1,w2,w3), and V1(w1,w2,w3) to V20(w1,w2,w3) the integration constants, which can be found in the Appendix (in Section A1.3). 3.2.2. The Gθ Functions. This set of FRFs corresponds to periodic modulation of the temperature of the adsorbent bed. The expressions for the first-, second-, and third-order FRFs are

first-order FRF: G1,θ(u1) ) H1(u1)eλ1(u1) + H2(u1)eλ2(u1) + h

(41)

The integration constants H1(u1), H2(u1), and h are given in Section A2.1 of the Appendix.


second-order FRF: G2,θθ(u1,u2) ) S1(u1,u2)e

+ S2(u1,u2)e

(λ1(u1)+λ1(u2))

+ s2(u1,u2)e

(λ2(u1)+λ1(u2))

+ s4(u1,u2)e

s1(u1,u2)e s3(u1,u2)e

µ1(u1,u2)

µ2(u1,u2)

+

(λ1(u1)+λ2(u2)) (λ2(u2)+λ2(u1))

+

third-order FRFs:

+

G3,ccθ(w1,w2,u1) ) Z1(w1,w2,u1)eχ1(w1,w2,u1) +

s5(u1,u2)eλ1(u1) + s6(u1,u2)eλ2(u2) + s7(u1,u2)eλ1(u2) + s8(u1,u2)e

λ2(u1)

The integration constants K1(w1,u1), K2(w1,u1) and k1(w1,u1) to k6(w1,u1) can be found in Section A3.1 of the Appendix.

+ s9 (42)

The integration constants S1(u1,u2), S2(u1,u2), and s1(u1,u2) to s9(u1,u2) can again be found in Section A2.2 of the Appendix.

Z2(w1,w2,u1)eχ2(w1,w2,u1) + z1(w1,w2,u1)e(µ1(w1,u1)+λ1(w2)) + z2(w1,w2,u1)e(µ2(w1,u1)+λ1(w2)) + z3(w1,w2,u1)e(µ1(w1,u1)+λ2(w2)) + z4(w1,w2,u1)e(µ2(w1,u1)+λ2(w2)) + z5(w1,w2,u1)e(µ1(w2,u1)+λ1(w1)) + z6(w1,w2,u1)e(µ2(w2,u1)+λ1(w1)) + z7(w1,w2,u1)e(µ1(w2,u1)+λ2(w1)) + z8(w1,w2,u1)e(µ2(w2,u1)+λ2(w1)) + z9(w1,w2,u1)e(µ1(w1,w2)+λ1(u1)) +

third-order FRF: G3,θθθ(u1,u2,u3) ) P1(u1,u2,u3)eχ1(u1,u2,u3) + P2(u1,u2,u3)eχ2(u1,u2,u3) + p1(u1,u2,u3)e(µ1(u1,u2)+λ1(u3)) + p2(u1,u2,u3)e(µ2(u1,u2)+λ1(u3)) + p3(u1,u2,u3)e(µ1(u1,u2)+λ2(u3)) + p4(u1,u2,u3)e(µ2(u1,u2)+λ2(u3)) + p5(u1,u2,u3)e(µ1(u1,u3)+λ1(u2)) + p6(u1,u2,u3)e(µ2(u1,u3)+λ1(u2)) + p7(u1,u2,u3)e(µ1(u1,u3)+λ2(u2)) + p8(u1,u2,u3)e(µ2(u1,u3)+λ2(u2)) + p9(u1,u2,u3)e(µ1(u2,u3)+λ1(u1)) + p10(u1,u2,u3)e(µ2(u2,u3)+λ1(u1)) + p11(u1,u2,u3)e(µ1(u2,u3)+λ2(u1)) + p12(u1,u2,u3)e(µ2(u2,u3)+λ2(u1)) + p13(u1,u2,u3)eµ1(u1,u2) + p14(u1,u2,u3)eµ2(u1,u2) + p15(u1,u2,u3)eµ1(u1,u3) + p16(u1,u2,u3)eµ2(u1,u3) + p17(u1,u2,u3)eµ1(u2,u3) + p18(u1,u2,u3)eµ2(u2,u3) + p19(u1,u2,u3)e(λ1(u1)+λ1(u2)+λ1(u3)) + p20(u1,u2,u3)e(λ1(u1)+λ2(u2)+λ1(u3)) + p21(u1,u2,u3)e(λ1(u2)+λ2(u1)+λ1(u3)) + p22(u1,u2,u3)e(λ1(u3)+λ2(u1)+λ2(u2)) + p23(u1,u2,u3)e(λ1(u1)+λ1(u2)+λ2(u3)) + p24(u1,u2,u3)e(λ1(u1)+λ2(u2)+λ2(u3)) + p25(u1,u2,u3)e(λ1(u2)+λ2(u1)+λ2(u3)) + p26(u1,u2,u3)e(λ2(u1)+λ2(u2)+λ2(u3)) + p27(u1,u2,u3)e(λ1(u1)+λ1(u3)) + p28(u1,u2,u3)e(λ1(u1)+λ1(u2)) + p29(u1,u2,u3)e(λ1(u2)+λ1(u3)) + p30(u1,u2,u3)e(λ2(u2)+λ1(u3)) + p31(u1,u2,u3)e(λ2(u1)+λ1(u3)) + p32(u1,u2,u3)e(λ1(u1)+λ2(u3)) + p33(u1,u2,u3)e(λ2(u2)+λ2(u3)) + p34(u1,u2,u3)e(λ1(u2)+λ2(u3)) + p35(u1,u2,u3)e(λ2(u1)+λ2(u3)) + p36(u1,u2,u3)e(λ1(u1)+λ2(u2)) + p37(u1,u2,u3)e(λ2(u1)+λ1(u2)) + p38(u1,u2,u3)e(λ2(u1)+λ2(u2)) + p39(u1,u2,u3)e(λ1(u1)) + p40(u1,u2,u3)e(λ1(u2)) + p41(u1,u2,u3)e(λ1(u3)) + p42(u1,u2,u3)e(λ2(u1)) + p43(u1,u2,u3)e(λ2(u2)) + p44(u1,u2,u3)e(λ2(u3)) + p45(u1,u2,u3) (43) The integration constants P1(u1,u2,u3), P2(u1,u2,u3), and p1(u1,u2,u3) to p45(u1,u2,u3) are also given in Section A2.3 of the Appendix. 3.2.3. The Gcθ Functions. This set contains the crossfunctions correlating the outlet concentration to both inputs (the inlet concentration and the adsorbent bed temperature). One second-order and two third-order functions can be defined:

second-order FRF: G2,cθ(w1,u1) ) K1(w1,u1)eµ1(w1,u1) + K2(w1,u1)eµ2(w1,u1) + k1(w1,u1)e(λ1(w1)+λ1(u1)) + k2(w1,u1)e(λ1(w1)+λ2(u1)) + k3(w1,u1)e(λ2(w1)+λ1(u1)) + k4(w1,u1)e(λ2(w1)+λ2(u1)) + k5(w1,u1)eλ1(w1) + k6(w1,u1)eλ2(w1) (44)

z10(w1,w2,u1)e(µ2(w1,w2)+λ1(u1)) + z11(w1,w2,u1)e(µ1(w1,w2)+λ2(u1)) + z12(w1,w2,u1)e(µ2(w1,w2)+λ2(u1)) + z13(w1,w2,u1)eµ1(w1,w2) + z14(w1,w2,u1)eµ2(w1,w2) + z15(w1,w2,u1)e(λ1(w1)+λ1(w2)+λ1(u1)) + z16(w1,w2,u1)e(λ1(w1)+λ1(w2)+λ2(u1)) + z17(w1,w2,u1)e(λ2(w1)+λ1(w2)+λ1(u1)) + z18(w1,w2,u1)e(λ1(w2)+λ2(w1)+λ2(u1)) + z19(w1,w2,u1)e(λ1(w1)+λ2(w2)+λ1(u1)) + z20(w1,w2,u1)e(λ1(w1)+λ2(w2)+λ2(u1)) + z21(w1,w2,u1)e(λ2(w1)+λ2(w2)+λ1(u1)) + z22(w1,w2,u1)e(λ2(w1)+λ2(w2)+λ2(u1)) + z23(w1,w2,u1)e(λ1(w1)+λ1(w2)) + z24(w1,w2,u1)e(λ1(w2)+λ2(w1)) + z25(w1,w2,u1)e(λ1(w1)+λ2(w2)) + z26(w1,w2,u1)e(λ2(w1)+λ2(w2)) (45) G3,cθθ(w1,u1,u2) ) R1(w1,u1,u2)eχ1(w1,u1,u2) + R2(w1,u1,u2)eχ2(w1,u1,u2) + r1(w1,u1,u2)e(µ1(u1,u2)+λ1(w1)) + r2(w1,u1,u2)e(µ2(u1,u2)+λ1(w1)) + r3(w1,u1,u2)e(µ1(u1,u2)+λ2(w1)) + r4(w1,u1,u2)e(µ2(u1,u2)+λ2(w1)) + r5(w1,u1,u2)e(µ1(w1,u2)+λ1(u1)) + r6(w1,u1,u2)e(µ2(w1,u2)+λ1(u1)) + r7(w1,u1,u2)e(µ1(w1,u2)+λ2(u1)) + r8(w1,u1,u2)e(µ2(w1,u2)+λ2(u1)) + r9(w1,u1,u2)e(µ1(w1,u1)+λ1(u2)) + r10(w1,u1,u2)e(µ2(w1,u1)+λ1(u2)) + r11(w1,u1,u2)e(µ1(w1,u1)+λ2(u2)) + r12(w1,u1,u2)e(µ2(w1,u1)+λ2(u2)) + r13(w1,u1,u2)eµ1(w1,u2) + r14(w1,u1,u2)eµ2(w1,u2) + r15(w1,u1,u2)eµ1(w1,u1) + r16(w1,u1,u2)eµ2(w1,u1) + r17(w1,u1,u2)e(λ1(w1)+λ1(u1)+λ1(u2)) + r18(w1,u1,u2)e(λ1(w1)+λ1(u1)+λ2(u2)) + r19(w1,u1,u2)e(λ2(w1)+λ1(u1)+λ1(u2)) + r20(w1,u1,u2)e(λ1(w1)+λ2(u1)+λ2(u2)) + r21(w,u1,u2)e(λ2(w1)+λ1(u1)+λ1(u2)) + r22(w1,u1,u2)e(λ2(w1)+λ2(u2)+λ1(u1)) + r23(w1,u1,u2)e(λ2(w1)+λ2(u1)+λ1(u2)) + (λ2(w1)+λ2(u1)+λ2(u2)) + r25(w1,u1,u2)e(λ1(w1)+λ1(u1)) + r24(w1,u1,u2)e r26(w1,u1,u2)e(λ1(w1)+λ2(u2)) + r27(w1,u1,u2)e(λ1(w1)+λ1(u2)) + r28(w,u1,u2)e(λ2(w)+λ2(u1)) + r29(w,u1,u2)e(λ2(w)+λ1(u1)) + r30(w,u1,u2)e(λ2(w)+λ2(u2)) + r31(w,u1,u2)e(λ2(w)+λ1(u2)) + r32(w1,u1,u2)e(λ2(w1)+λ2(u1)) + r33(w1,u1,u2)e(λ1(w1)) + r34(w1,u1,u2)e(λ2(w1)) (46) The integration constants Z1(w1,w2,u1), Z2(w1,w2,u1), and z1(w1,w2,u1) to z26(w1,w2,u1), as well as R1(w1,u1,u2), R2(w1,u1,u2)

272


Figure 3. Simulated first-, second-, and third-order Gc functions: G1,c(ω), G2,cc(ω,ω), and G3,ccc(ω,ω,ω).

Figure 4. Simulated first-, second-, and third-order Gθ functions: G1,θ(ω), G2,θθ(ω,ω), and G3,θθθ(ω,ω,ω).

and r1(w1,u1,u2) to r34(w1,u1,u2), can be found in Section A3.2 of the Appendix. Note that all frequencies used in the FRFs expressions (w1 through w3 and u1 through u3) are nondimensional. 3.2.4. Simulation of Adsorption Column Frequency Response Functions. To illustrate the FRFs of the adsorption column, given in the previous sections, the derived expressions were used for their simulation. The simulation results presented here were obtained for a laboratory-scale column of the following dimensions and characteristics: length, L ) 10 cm; velocity, V ) 0.0737 cm/s; ) 0.602; NTP ) 1000; and for a Langmuir adsorption isotherm described by

Q)

Q0bC 1 + bC

(Tk)

b ) b0 exp

with the coefficients Q0 ) 0.0885 mol/dm3, b0 ) 0.00765 dm3/ mol, and k ) 2647 K-1. These parameter values correspond to the following system: microcrystalline cellulose triacetate as the adsorbent, Tro¨ger’s base (-) enanthiomer as the adsorbate, and ethanol as the carrier fluid.27 The isotherm derivatives in the expressions for the FRFs were calculated using eqs 37 for the following steady-state concentration and temperature: Cs ) 0.005 mol/dm3 and Ts ) 315 K. The three sets of the simulated HFRFs, up to the third order, are presented in Figures 3 (the Gc functions), 4 (the Gθ functions), and 5 (the Gcθ functions). All these functions correspond to the case u1 ) u2 ) u3 ) w1 ) w2 ) w3 ) ω, i.e., they correspond to the dominant terms of the first, second, and third harmonic of the nondimensional outlet concentration co, for single harmonic changes of the nondimensional inlet concentration into the column ci or/and nondimensional temperature of the column θ, of the same frequency ω. All FRFs are presented in standard Bode´ plot form (amplitude functions versus frequency in log-log diagrams, and phase functions versus frequency in semilogarithmic diagrams). Without discussing the characteristics of all the presented FRFs in detail, we will make a note of only a couple of the most significant issues: All FRFs except G1,c(ω) approach a value of zero at low frequencies, whereas G1,c(ω) approaches a value of 1. As a

Figure 5. Simulated second- and third-order Gcθ functions: G2,cθ(ω,ω), G3,ccθ(ω,ω,ω), and G3,cθθ(ω,ω,ω).

consequence, for very slow modulations of the inputs, the outlet concentration from the column follows the change of the inlet concentration, regardless of the temperature change. Certain oscilatory behavior of the amplitudes of the Gθ and Gcθ functions can be observed, with a great number of resonant frequencies and a great number of frequencies for which the amplitudes approach a value of zero. This behavior is not unexpected, because these FRFs correspond to periodic changes in the temperature, which is introduced to the entire adsorbent bed, i.e., in a distributed manner. This behavior has already been reported for shell-and-tube heat exchangers with a sinusoidal change of the heating fluid temperature.5 3.3. Fast Approximate Calculation of the Outlet Concentration in the Periodic Quasi-Steady State. Using the HFRFs derived in the previous section, the outlet concentration from the adsorption column in the periodic quasi-steady state can be calculated approximately. As mentioned in Section 2.2, contrary to the exact response of a nonlinear system, which contains a DC component and an indefinite number of harmonics of different frequencies, only a finite number of harmonics appear in the approximate solution.


Figure 6. Inlet concentration and adsorbent temperature input changes (under conditions of Cs ) 0.005 mol/dm3, Ts ) 315 K, A ) 25%, B ) 3.2%).

Figure 7. Outlet concentrations and their Fourier spectra for 25% cosinusoidal change of the inlet concentration, for Case 1 and a period of (a) 162.5 s and (b) 628 s.

Let us consider a periodic behavior of an adsorption column with arbitrary shapes of the input periodic changes of the inlet concentration or/and the temperature of the adsorbent bed. By applying the Fourier analysis and taking into consideration only the first K harmonics, these periodic inputs can be approximated by finite sums of simple harmonic functions: K

ci,app(t) )

∑|Ak| cos(kωt + arg(Ak))

(47)

In eqs 47 and 48, the basic frequency of both periodic inputs was chosen to be the same. Ak and Bk are the kth members of the Fourier transforms of the two inputs. Their complex magnitudes define the amplitudes, and their phase angles define the phases, of the kth harmonics of the inputs. If the nonlinear model of the adsorption column is approximated by a set of N HFRFs, the approximate outlet concentration is obtained as a sum of a DC component and N × K harmonics:

k)1 K

θapp(t) )

∑|Bk| cos(kωt + arg(Bk)) k)1

K×N

(48)

co,app(t) )

∑ |Yn| cos(nωt + arg(Yn))

n)0

(49)

274


Figure 8. Outlet concentrations and their Fourier spectra for 3.2% cosinusoidal change of the adsorbent temperature, for Case 2 and a period of (a) 162.5 s and (b) 628 s.

The approximate expressions for the nth harmonic, of frequency nω, in its complex form, are as follows: (1) For Case 1 (periodic change of the inlet concentration only): K

Yn(nω) )

K

1

K

1

G2,cc(k1ω,k2ω)Ak1Ak2 + ‚‚‚ +

K

K

∑ ∑

2N-1 k1)-K k2)-K

K

‚‚‚

∑

kN)-K |k1+k2+...kN|)n

GN,cN (k1ω,k2ω,...,kNω)Ak1Ak2‚‚‚AkN (50) (2) For Case 2 (periodic change of the adsorbent temperature only): K

Yn(nω) )

∑

k)-K |k|)n

G1,θ(kω)Bk +

1

K

K

∑ ∑

2 k1)-K k2)-K |k1+k2|)n

G2,θθ(k1ω,k2ω)Bk1Bk2 + ‚‚‚ +

K

1

K

K

∑ ∑ ‚‚‚ k ∑ N-1 k )-K k )-Kx )-K

2

1

Yn(nω) )

K

1

K

∑ (G1,c(kω)Ak + G1,θ(kω)Bk) + 2 k ∑ ∑ )-K k )-K 1

2

|k1+k2|)n

(G2,cc(k1ω,k2ω)Ak2Ak2 + G2,cθ(k1ω,k2ω)Ak1Bk2 +

2

|k1+k2|)n

1

K

k)-K |k|)n

∑ G1,c(kω)Ak + 2 k ∑ ∑ k)-K )-K k )-K |k|)n

(3) For Case 3 (periodic changes of both the inlet concentration and the adsorbent temperature):

2

N

|k1+k2+...kN|)n

GN,θN(k1ω,k2ω,...,kNω)Bk1Bk2‚‚‚BkN (51)

G2,θθ(k1ω,k2ω)Bk1Bk2) + ‚‚‚ +

1

N

K

K

K

2

∑ ∑ ∑ ‚‚‚ k ∑ N-1 m)0 k )-K k )-K )-K 1

2

N

|k1+k2+...kN|)n

GN,cmθN-m(k1ω,k2ω,...,kNω)Ak2‚‚‚AkmBkm+1‚‚‚BkN (52) 3.3.1. Example 1. Adsorption Column Response to Single Harmonic Inputs. To illustrate the method for fast approximate calculation of periodic quasi-steady states of an adsorption column, presented in the previous text, we will first apply it to the simplest case of single harmonic (cosinusoidal) input changes. For this case, K ) 1 in eqs 47-52. The HFRFs, up to the third order, which have been derived in the Appendix and listed in Section 3.2, are used; i.e., in eqs 47-52, N ) 3. A sample of the results obtained for the test system defined in Section 3.2.4 and for the input functions presented in Figure 6 (25% nondimensional input concentration change of amplitude and 3.2% nondimensional temperature change around the steady state defined by Cs ) 0.005 mol/dm3 and Ts ) 315 K) is given


Figure 9. Outlet concentrations and their Fourier spectra for simultaneous, out-of-phase cosinusoidal changes of the inlet concentration (25%) and adsorbent temperature (3.2%), for Case 3 and a period of (a) 162.5 s and (b) 628 s.

in Figures 7, 8, and 9, for Cases 1, 2, and 3, respectively. For Case 3, the changes of the inlet concentration and temperature are out of phase, which corresponds to real cyclic adsorption processes: adsorption at low and desorption at high temperatures. Note that a 3.2% amplitude of temperature change corresponds to a difference of ∼20 K between the maximal and minimal temperature. This temperature change results in approximately the same change of the concentration in the solid phase as 25% change of the concentration in the fluid phase. For each case, we give the results obtained for two different periods: 162.5 and 628 s, corresponding to nondimensional frequencies of 5.2464 and 1.3575, respectively. These frequencies were chosen in such a way to correspond to different behavior of the FRFs: high nonlinearity with distinctive periodic behavior of the Gθ and Gcθ functions (162.5 s) and moderate nonlinearity (628 s). For comparison, along with the approximate outlet concentration calculated based on the first three HFRFs, we also give the numerical solution of the model equations described by eqs 29-33, which is assumed to be exact. The numerical solution was obtained using the fast converging algorithm of SeidelMorgenstern et al.28 A segment of quasi-steady-state portion of the numerical solution was used for comparison with the approximate solution in Figures 7-9. Along with the outlet concentrations in time domain, we give the Fourier spectra of the output, i.e., the amplitudes and phases of the Fourier transforms of the approximate and the exact solutions. The amplitudes correspond to the nondimensional outlet concentrations.

Figures 7-9 show that, even with only three HFRFs, very good agreement of the approximate and the exact solution for the outlet concentration is obtained. This agreement is especially good for Case 1, when only the inlet concentration is modulated. The agreement is good in regard to both the time and the frequency domain. 3.3.2. Example 2. Adsorption Column Response to Rectangular Inputs. In practice, the periodic adsorption/desorption is not performed with sinusoidal, but rather with rectangular inlet concentration and temperature changes. Our method for fast approximate calculation of the outlet concentration was tested for such a case as well. In Figure 10, we show the input changes of the rectangular shape, together with their Fourier spectra. Again, for Case 3, the inlet concentration and the temperature waves are out of phase. In practice, the desorption/ hot portion of the cycle is considerably shorter than the adsorption/cold portion. Such is also the situation in our simulations. The ratio of the durations of the adsorption and desorption half-cycles is 3:2. It is interesting to note that, for this ratio, every fifth harmonic in the inputs vanishes. The inputs shown in Figure 10 correspond to a steady-state concentration and temperature of 0.005 mol/dm3 and 315 K, respectively, and for inlet concentration and temperature amplitudes of 25% and 3.2%, respectively. The approximate concentrations in the periodic quasi-steady state corresponding to the inputs shown in Figure 10, and for Cases 1, 2, and 3 are shown in Figures 11, 12, and 13, respectively. The calculations were performed for the same

276


Figure 10. Rectangular input changes and their Fourier spectra: (a) 25% change of the inlet concentration, and (b) 3.2% change of the adsorbent temperature. Table 1. Relative Errors of the Average Concentrations during the Adsorption and Desorption Half-Cycle, for Cosinusoidal Input Changes Case 1 A ) 12.5%

A ) 25.0%

A ) 50.0%

A ) 100.0%

period (s)

N

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

162.5 162.5 162.5 628 628 628

1 2 3 1 2 3

-0.499 -0.355 -0.148 -0.072 -0.047 -0.029

0.499 0.354 0.151 0.082 0.060 0.033

1.474 0.936 0.506 -0.136 -0.113 -0.093

-1.471 -0.946 -0.637 0.157 0.095 0.087 Case 2

-2.299 -2.052 0.948 0.314 0.220 0.127

1.799 1.690 -0.960 -0.323 -0.222 -0.121

16.907 -8.502 -6.015 0.654 0.306 0.281

-14.609 6.372 4.950 -0.792 -0.312 -0.290

period (s)

N

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

162.5 162.5 162.5 628 628 628

1 2 3 1 2 3

-0.902 -0.811 -0.534 0.129 0.110 0.073

0.891 0.788 0.515 -0.126 -0.975 -0.069

-1.312 -1.208 -0.904 0.310 0.208 0.157

1.302 1.197 0.805 -0.304 -0.205 -0.153 Case 3

3.861 2.749 2.166 -1.609 -1.552 -0.725

-3.912 -3.125 -2.669 1.679 1.591 0.756

15.208 -10.594 -9.727 7.733 5.788 2.223

-14.257 9.469 8.348 -7.637 -5.583 -2.193

B ) 1.6%

B ) 3.2%

B ) 8.0%

B ) 15.0%

A ) 12.5%, B ) 1.6%

A ) 25.0%, B ) 3.2%

A ) 50.0%, B ) 8.0%

period (s)

N

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

162.5 162.5 162.5 628 628 628

1 2 3 1 2 3

-1.381 -1.261 -0.950 0.259 0.202 0.188

1.352 1.299 0.971 -0.212 -0.198 -0.183

-2.845 -1.999 -1.950 0.857 0.831 0.717

2.499 1.846 1.700 -0.858 -0.840 -0.722

-20.662 -11.219 -10.989 6.585 5.432 3.104

17.991 9.578 9.186 -6.390 -5.293 -3.974

-36.341 -25.719 -20.275 14.729 12.328 7.277

30.681 24.275 18.722 -12.107 -8.384 -5.620

periods as for the single harmonic inputs (162.5 and 628 s). In these calculations, the HFRFs up to the third order were used again (N ) 3), and four harmonics of the inputs were taken into account (K ) 4). The numerical (exact) solution is also

A ) 100.0%, B ) 15.0%

given, as well as the frequency spectra corresponding to the approximate and exact solutions. 3.5. Analysisscomparison with Numerical Solution. The results presented in the previous section show relatively good


Figure 11. Outlet concentrations and their Fourier spectra for 25% rectangular change of the inlet concentration, for Case 1 and a period of (a) 162.5 s and (b) 628 s.

agreement between the approximate (based on the first three FRFs) and exact (numerical) solution of the adsorption column quasi-steady state response, for all three cases under consideration. One can observe that better agreement was obtained for Case 1, than for the other two cases. Also, the agreement is better for cosinusoidal than for rectangular inputs, which was expected, as approximative input functions with only four harmonics were used in the second case. Regarding the input period, better agreement was obtained for 628 s than for 162.5 s. This is also an expected result, because a period of 162.5 s corresponds to very distinctive nonlinear behavior of the adsorption column, especially with regard to the change of the column temperature. In real applications, the exact time profiles of the outlet concentration are not of great importance; however, their average values during the adsorption and desorption half-cycles are. For that reason, these average values were calculated using the approximate solution and the exact solution, and the quality of the approximate solution was estimated based on the relative errors defined in the following way:

δA(%) ) δD(%) )

〈Co,app〉A - 〈Co,ex〉A 〈Co,ex〉A 〈Co,app〉D - 〈Co,ex〉D 〈Co,ex〉D

× 100

(53)

× 100

(54)

In these equations, the use of angled brackets (〈〉) denotes the

average values, the subscript app corresponds to the approximate solution (based on a finite number of HFRFs), ex refers to the exact (numerical) solution, A refers to the adsorption half-cycle, and D refers to the desorption half-cycle. The average values of the outlet concentration were calculated by numerical integration. Table 1 gives the relative errors, defined by eqs 53 and 54, obtained for cosinusoidal input changes, for the test system defined in Section 3.2.4 and for Cs ) 0.005 mol/dm3 and Ts ) 315 K. The relative errors obtained for Cases 1, 2, and 3, for four different combinations of input amplitudes, and for two periods, are listed. The concentration and temperature amplitudes were chosen in such a way to cause approximately the same change of solid concentration: a 12.5% change of concentration causes approximately the same change of solid concentration as a 1.6% change of temperature, a 25% change of concentration corresponds to an ∼3.2% change of temperature, a 50% change of concentration relates to an ∼8% change of temperature, and a 100% change of concentration corresponds to an ∼15% change of temperature. High input amplitudes were chosen for analysis, to check the applicability of our method for fast approximate calculation of the quasi-steady states for real applications. To analyze the influence of the number of HFRFs used for calculation of the approximate solution, in Table 1, the relative errors obtained with N ) 1, 2, or 3 FRFs are given for each case. In Table 2, we give the relative errors obtained for the case of rectangular input waves. The input amplitudes and periods

278


Figure 12. Outlet concentrations and their Fourier spectra for 3.2% rectangular change of the adsorbent bed temperature, for Case 2 and a period of (a) 162.5 s and (b) 628 s.

chosen for analysis were the same as those used in Table 1. All presented results were obtained for N ) 3 (three FRFs were used for the calculation of the approximate solution); however, two different values of the number of input harmonics that were taken into considerationsK ) 4 and K ) 9sare used. By inspection of the relative errors given in Tables 1 and 2, the following can be concluded: (1) The relative errors decrease as more FR functions are included (increasing N); (2) The relative errors decrease as more input harmonics are used (increasing K), in the case of rectangular input waves; (3) The relative errors increase as the input amplitude(s) increases; (4) The relative errors are larger for the period of 162.5 s than for the period of 628 s; (5) The smallest relative errors are obtained for Case 1 and the largest are obtained for Case 3; (6) Somewhat smaller relative errors were obtained with cosinusoidal input waves than with rectangular input waves. All of these results were expected. 4. Conclusions The new method for fast estimation of periodic quasi-steady states of cyclic processes, presented in this paper, was developed with the idea to substitute long numerical simulations. The method is based on Volterra series approximation and the higher-order frequency response functions (HFRFs), which are derived starting from a nonlinear model of the process.

The complete procedure used for application of this method consist of four consecutive steps: postulating of an appropriate mathematical model, derivation of the HFRFs, approximation of the inputs by finite sums, and calculation of the approximate output. It is important to note that the most time-consuming step of this procedure, which is the derivation of the frequency response functions (FRFs), needs to be performed only once. Although the mathematical expressions used in the other steps appear very cumbersome, they actually reduce to simple algebra. In the proposed method, the nonlinear system model is approximated by a sequence of a finite number of FRFs and the input modulation functions are approximated by sums of a finite number of harmonics. As a consequence, the proposed method is used for approximate calculation of periodic change of the system output in its quasi-steady state. The quality of the approximate solution obtained by the proposed method, i.e., its closeness to the exact solution, is directly dependent on the number of HFRFs and the number of input harmonics used in the calculation. Nevertheless, the method enables direct calculation of the quasi-steady state, by using simple algebraic operations. On the other hand, the commonly used numerical solutions demand calculations of a large number of cycles to obtain negligible transient response and quasi-steady state and involve calculations that are much more complicated and time-consuming. Application of the proposed method was illustrated using a relatively simple case study: the cyclic operation of an adsorption column with periodic modulation of the inlet


Figure 13. Outlet concentrations and their Fourier spectra for simultaneous out-of-phase rectangular changes of the inlet concentration (25%) and adsorbent bed temperature (3.2%), for Case 3 and a period of (a) 162.5 s and (b) 628 s. Table 2. Relative Errors of the Average Concentrations during the Adsorption and Desorption Half-Cycle, for Rectangular Input Changes Case 1 A ) 12.5%

A ) 25.0%

A ) 50.0%

A )100.0%

period (s)

K

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

162.5 162.5 628 628

4 9 4 9

0.417 0.405 0.069 0.013

-0.272 -0.233 -0.081 -0.052

0.660 0.593 0.202 0.030

-0.389 -0.322 -0.205 -0.066

2.796 1.498 0.278 0.143

-0.809 -0.365 -0.175 -0.053

6.094 4.286 0.700 0.242

-8.910 -5.622 -0.195 -0.176

Case 2 B ) 1.6%

B ) 3.2%

B ) 8.0%

B )15.0%

period (s)

K

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

162.5 162.5 628 628

4 9 4 9

-0.863 -0.504 -0.403 -0.210

0.607 0.352 0.366 0.193

-1.752 -1.083 -0.532 -0.388

1.316 0.813 0.636 0.445

-5.305 -3.770 -1.943 -1.691

5.413 3.493 1.134 0.963

-13.862 -11.390 -6.709 -4.546

20.527 19.001 8.728 6.188

Case 3 A ) 12.5%, B ) 1.6%

A ) 25.0%, B ) 3.2%

period (s)

K

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

δA (%)

δD (%)

162.5 162.5 628 628

4 9 4 9

-1.155 -0.751 -0.737 -0.537

0.687 0.439 0.784 0.577

-2.873 -1.904 -1.405 -0.919

3.179 2.076 0.864 0.740

14.405 12.316 5.344 4.108

-9.648 -8.001 -6.068 -4.893

-24.114 -21.145 9.163 7.993

21.690 19.827 -8.635 -7.028

concentration or/and temperature of the entire column. The approximate quasi-steady state periodic changes of the outlet concentrations, calculated using only the first-, second-, and third-order FRFs were calculated for cosinusoidal and rectangular changes of only one and both inputs. The quality of the results was estimated based on their agreement with the outlet concentrations calculated by numerical integration of the model

A ) 50.0%, B ) 8.0%

A ) 100.0%, B ) 15.0%

equations. The agreement was generally good, even for very large input amplitudes: 100% change of the inlet concentration, and 15% change of the column temperature (∼80 K between the maximal and minimal temperature during the cycle). Based on the results obtained on this case study, it can be concluded that the method has good potential for estimation of periodic behavior of different cyclic processes.

280


higher-order frequency response functions (HFRFs), which will be denoted as Fc, Fθ, and Fcθ functions, and which correlate the concentration in the fluid phase c at any position x in the column to the same inputs. For x ) 1, the F functions become identical to the G functions, which are the goals of the derivation procedure. After setting up the model equations in the nondimensional form (eqs 29-33 in Section 3.1), the HFRFs derivation is performed in five steps, which are defined in Section 3.2: Step 1. Defining the inputs. To derive the HFRFs, up to the third order, it is most convenient to define the inputs in the following way:

Finally, we must note that the quality of the results obtained with our method is dependent not only on the number of HFRFs and input harmonics used in the approximation, but primarily on the realism of the nonlinear model used for derivation of the HFRFs. Acknowledgment This work was supported by the Serbian Ministry of Science and Environmental Protection in the frame of Project No. 1700. The authors would like to thank Prof. Andreas Seidel-Morgenstern (Otto-von-Guericke University, Magdeburg, Germany) for useful discussions and for providing us with the computer code for the numerical solution.

3

Appendix: Derivation of the Adsorption Column Higher-Order Frequency Response Functions, up to the Third Order The Gc, Gθ, and Gcθ HFRFs that have been defined and used in the main body of this manucsript relate the change of the outlet concentration from the column to the modulation of the inlet concentration or/and the adsorbent temperature. To derive these functions, it is convenient to define three sets of auxiliary 3

∑

c(x,t) )

3

AiF1,c(x,wi)ejwit +

i)1 3

∑

3

i)1

3

3

∑∑A B F

i j 2,cθ(x,wi,uj)e

j(wi+uj)t

i)1 j)1

+

(A-1)

3

Bkeju t ∑ k)1 k

(A-2)

Step 2. Representing the concentration in the fluid phase at position x, in the Volterra series form: 3

∑∑∑A A A F

i j k 3,ccc(x,wi,wj,wk)e

3

∑∑B B F

AiAjF2,cc(x,wi,wj)ej(wi+wj)t +

i)1 j)1 3

3

i

θ(t) )

3

∑∑

BiF1,θ(x,ui)ejuit +

Aiejw t ∑ i)1

ci(t) )

j(ui+uj)t i j 2,θθ(x,ui,uj)e

i)1 j)1 3 (wi+wj+wk)t

+

i)1 j)1 k)1 3 3 3

3

+

3

∑∑∑B B B F

i j k 3,θθθ(x,ui,uj,uk)e

j(ui+uj+uk)t

+

i)1 j)1 k)1 3 3 3

∑∑∑ A A B F

j(wi+wj+uk)t i j k 3,ccθ(x,wi,wj,uk)e

+

i)1 j)1 k)1

∑∑∑ A B B F

j(wi+uj+uk)t i j k 3,cθθ(x,wi,uj,uk)e

+ ‚‚‚ (A-3)

i)1 j)1 k)1

Step 3. Expressing the concentration in the solid phase at position x, using the assumption of local equilibrium and the equilibrium relation (eq 32): 3

∑

q(x,t) ) a˜c

3

AiF1,c(x,wi)ejwit + a˜c

i)1

∑

3

∑∑

i)1

3

3

i j 2,cθ(x,wi,uj)e

j(wi+uj)t

3

3

3

i j k 3,ccθ(x,wi,wj,uk)e

(wi+wj+uk)t

3

∑∑

3

AiAjF1,c(x,wi)F1,c(x,wj)ej(wi+wj)t + b˜ cc

i)1 j)1

+ a˜c

3

3

3

3

3

3

3

3

j(wi+uj)t i j 1,c(x,wi)e

+ b˜ cθ

i)1 j)1

3

3

3

3

∑∑A B F

i j 1,c(x,wi)F1,θ(x,uj)e

(ui+uj+uk)t

+ b˜ cθ

i)1 j)1 k)1

i j k 2,θθ(x,ui,uj)F1,θ(x,uk)e

(ui+uj+uk)t

+

3

3

∑∑∑

AiBjBkF1,c(x,wi)F2,θθ(x,uj,uk)e(wi+uj+uk)t + ‚‚‚ + b˜ θθ

3

3

i j 1,θ(x,ui)e

j(ui+uj)t

+ b˜ cθ

3

3

i j k 2,cc(x,wi,wj)e 3

i j k 2,cθ(x,wi,uj)e

(wi+uj+uk)t

3

+ ‚‚‚ + c˜ccc 3

(wi+wj+uk)t

+

i j k 1,c(x,wi)F1,θ(x,uj)F1,θ(x,uk)e

3

3

∑∑∑A A A F

i j k 1,c(x,wi)F1,c(x,wj)

3

∑∑∑ A A B F

i j k 1,c(x,wi)F1,c(x,wj)F1,θ(x,uk)e

i)1 j)1 k)1 3 3

∑∑∑A B B F

(wi+uj+uk)t

+ ... + c˜θθθ

(wi+wj+uk)t

+

3

∑∑∑B B B e

j(ui+uj+uk)t

i j k

+

i)1 j)1 k)1

3

3

j(wi+wj+uk)t

+

i)1 j)1 k)1

BiBjBkF1,θ(x,ui)F1,θ(x,uj)F1,θ(x,uk)e(ui+uj+uk)t + c˜ccc

i j k 1,c(x,wi)F1,c(x,wj)e

j(ui+uj)t

i j

3

3

∑∑∑A B B F

i)1 j)1 k)1

∑∑∑A A B F

3

∑∑∑ A A B F

3

∑∑∑

∑∑B B e

i)1 j)1 k)1

i)1 j)1 k)1 3 3 3

c˜ccc

+

i)1 j)1

i)1 j)1 k)1

3

j(wi+uj)t

3

3

∑∑B B F 3

i j k 2,θθ(x,ui,uj)e

F1,c(x,wk)e(wi+wj+wk)t + c˜ccc

+

3

∑∑∑ B B B F

i)1 j)1

∑∑∑B B B F

b˜ cθ

3

juit

i

i)1 j)1 k)1

3

3

∑B e

i)1 j)1

3

∑∑A B F

c˜ccθ

+ ‚‚‚ + a˜θ

i)1 j)1 k)1

i)1 j)1 k)1

3

(wi+uj+uk)t

3

3

AiAjBkF2,cc(x,wi,wj)F1,θ(x,uk)e(wi+wj+uk)t + b˜ cc b˜ cθ

+

i)1

BiBjF1,θ(x,ui)F1,θ(x,uj)ej(ui+uj)t + b˜ cc

AiAjAkF2,cc(x,wi,wj)F1,c(x,wk)e(wi+wj+wk)t + b˜ cc

3

∑∑∑

(ui+uj+uk)t

3

i j k 3,cθθ(x,wi,uj,uk)e

i)1 j)1 k)1

b˜ cc

i j k 3,θθθ(x,ui,uj,uk)e

3

∑∑∑ A B B F

3

∑∑∑

3

i)1 j)1 k)1

i)1 j)1

b˜ cc

3

+ a˜c

3

∑∑

+

∑∑∑B B B F

(wi+wj+wk)t

i)1 j)1 k)1

3

j(ui+uj)t

i)1 j)1

i j k 3,ccc(x,wi,wj,wk)e

i)1 j)1 k)1

b˜ cc

i j 2,θθ(x,ui,uj)e

i)1 j)1 k)1

∑∑∑A A B F

a˜c

3

∑∑B B F

∑∑∑A A A F

+ a˜c

i)1 j)1 3

3

AiAjF2,cc(x,wi,wj)ej(wi+wj)t + a˜c

i)1 j)1 3 3

3

∑∑A B F

a˜c

3

BiF1,θ(x,ui)ejuit + a˜c

+ ... + c˜ccθ

i)1 j)1 k)1

3

3

j(wi+uj+uk)t i j k 1,c(x,wi)F1,θ(x,uj)e

3

3

∑∑∑ i)1 j)1 k)1

3

i j k 1,θ

i)1 j)1 k)1

3

AiBjBkF1,c(x,wi)ej(wi+uj+uk)t + c˜cθθ

3

∑∑∑B B B F

+ ‚‚‚ + c˜ccθ

i)1 j)1 k)1

3

(x,ui)F1,θ(x,uj)ej(ui+uj+uk)t + ‚‚‚ + c˜cθθ

3

∑∑∑ A B B F

3

3

∑∑∑B B B F

i j k 1,θ(x,ui)e

i)1 j)1 k)1

j(ui+uj+uk)t

+ ... (A-4)


Step 4. Substitution of the expressions for ci, θ, c, and q into the model eq 29 and boundary conditions (eqs 30 and 31). The substitution itself is trivial, but the resulting equations are very cumbersome, so they are not listed here. Step 5. Applying the method of harmonic probing to the equations obtained in step 4 and solution of the resulting equations. This step will be shown for each derived function separately. A1. The Gc Functions. A1.1. The First-Order Gc Function. By collecting the terms containing A1ejw1t in the model equations obtained in Step 4, and equating them to 0, the following linear, homogeneous second-order ordinary differential equation (ODE) is obtained:

d2F1,c(x,w1) dx

2

- 2NTP

dF1,c(x,w1) dx 2NTPjw1F1,c(x,w1)(1 + a˜cf) ) 0 (A-5)

The solution is

F2,cc(x,w1,w2) ) D1(w1,w2)eµ1(w1,w2)x + D2(w1,w2)eµ2(w1,w2)x + d1(w1,w2)e(λ1(w1) + λ1(w2))x + d2(w1,w2)e(λ1(w1)+λ2(w2))x + d3(w1,w2)e(λ2(w1)+λ1(w2))x + d4(w1,w2)e(λ2(w1)+λ2(w2))x (A-15) where µ1,2, defined as

µ1,2(w1,w2) ) NTP ( xNTP2 + 2NTPj(1 + a˜cf)(w1 + w2) (A-16) are the characteristic values, d1(w1,w2) through d4(w1,w2) are the integration constants corresponding to the particular solution:

with the following boundary conditions:

x ) 0: F1,c(0,w1) ) 1 + x ) 1:

1 dF1,c(x,w1) | | (A-6) 2NTP dx |x)0

dF1,c(x,w1) | | )0 dx |x)1

(A-7)

The solution of eq A-5 is

F1,c(x,w1) ) U1(w1)eλ1(w1)x + U2(w1)eλ2(w1)x

(A-9)

and U1(w1) and U2(w1) are the integration constants obtained from the following set of equations resulting from the boundary conditions described by eqs A-6 and A-7:

U1(w1) + U2(w1) ) 1 +

(U1(w1)λ1(w1) + U2(w1)λ2(w1))/2 NTP (A-10)

U1(w1)λ1(w1)eλ1(w1) + U2(w1)λ2(w1)eλ2(w1) ) 0

d2 )

NTPj(w1 + w2) fb˜ ccU1,1U2,2 λ1,1λ2,2

d3 )

NTPj(w1 + w2) fb˜ ccU1,2U2,1 λ1,2λ2,1

d4 )

NTPj(w1 + w2) fb˜ ccU2,1U2,2 λ2,1λ2,2

and the integration constants D1(w1,w2) and D2(w1,w2), corresponding to the complementary solution, are obtained from the following set of algebraic equations:

D1 + D2 + d1 + d2 + d3 + d4 ) {[µ1D1 + µ2D2 + d1(λ1,1 + λ1,2) + d2(λ1,1 + λ2,2) + d3(λ2,1 + λ1,2) + d4(λ2,1 + λ2,2)]/2}/NTP (A-17) D1 µ1e µ1 + D2 µ2e µ2 + d1(λ1,1 + λ1,2)e(λ1,1+λ1,2) + d2(λ1,1 + λ2,2)e(λ1,1+λ2,2) + d3(λ2,1 + λ1,2)e(λ2,1+λ1,2) + d4(λ2,1 + λ2,2)e(λ2,1+λ2,2) ) 0 (A-18)

(A-11)

Because, by definition, G1,c(w1) ≡ F1,c(x ) 1, w1), by substituting x ) 1 into eq A-8, the final expression for the first-order Gc function (eq 38 in Section 3.2.1) is obtained. A1.2. The Second-Order Gc Function. This function is obtained by collecting the terms of the equations obtained in Step 4 that contain A1A2ej(w1+w2)t and then equating them to zero. In this way, the model equation is transformed to the following linear, nonhomogeneous, second-order ODE:

d2F2,cc(x,w1,w2)

dF2,cc(x,w1,w2) - 2NTP dx dx 2NTPj(w1 + w2)F2,cc(x,w1,w2)(1 + a˜c f) ) 2NTPj(w1 + w2) f b˜ ccF1,c(x,w1)F1,c(x,w2) (A-12) 2

and the boundary conditions become

The following shortened notations have been used to simplify the expressions: di ) di(w1,w2), Di ) Di(w1,w2), µi ) µi(w1,w2), Ui,j ) Ui(wj), and λi,j ) λi(wj). For x ) 1, eq A-15 reduces to eq 39 in Section 3.2.1, defining the function G2,cc(w1,w2). A1.3. The Third-Order Gc Function. This function is obtained by collecting the terms of the equations obtained in Step 4 that contain A1A2A3ej(w1+w2+w3)t and equating them to zero. The following linear, nonhomogeneous, second-order ODE is obtained:

d2F3,ccc(x,w1,w2,w3)

dF3,ccc(x,w1,w2,w3) - 2NTP dx dx 2NTP j(w1 + w2 + w3)F3,ccc(x,w1,w2,w3)(1 + a˜c f) ) 2

[32 b˜ (F cc

1 dF2,cc(x,w1,w2) | x ) 0: F2,cc(0,w1,w2) ) | (A-13) 2NTP dx |x)0 dF2,cc(x,w1,w2) | | )0 x ) 1: dx |x)1

NTPj(w1 + w2) fb˜ ccU1,1U1,2 λ1,1λ1,2

(A-8)

where λ1(w1)and λ2(w1) are its characteristic values,

λ1,2(w1) ) NTP ( xNTP2 + 2NTP(1 + a˜cf)jw1

d1 )

(A-14)

2,cc(x,w1,w2)F1,c(x,w3)

+ F2,cc(x,w1,w3)F1,c(x,w2) +

F2,c(x,w2,w3)F1,c(x,w1)) + c˜cccF1,c(x,w1)F1,c(x,w2) ×

]

F1,c(x,w3) 2NTP j(w1 + w2 + w3)f (A-19) with the following boundary conditions:

282


1 dF3,ccc(x,w1,w2,w3) | | 2NTP dx |x)0 (A-20) dF3,ccc(x,w1,w2,w3) | | )0 (A-21) x ) 1: dx |x)1

x ) 0: F3,ccc(0,w1,w2,w3) )

The solution is F3,ccc(x,w1,w2,w3) ) V1(w1,w2,w3)eχ1(w1,w2,w3)x + V2(w1,w2,w3)eχ2(w1,w2,w3)x + V1(w1,w2,w3)e(µ1(w1,w2)+λ1(w3))x + V2(w1,w2,w3)e(µ2(w1,w2)+λ1(w3))x + V3(w1,w2,w3)e(µ1(w1,w2)+λ2(w3))x + V4(w1,w2,w3)e(µ2(w1,w2)+λ2(w3))x + V5(w1,w2,w3)e(µ1(w1,w3)+λ1(w2))x + V6(w1,w2,w3)e(µ2(w1,w3)+λ1(w2))x + V7(w1,w2,w3)e(µ1(w1,w3)+λ2(w2))x + V8(w1,w2,w3)e(µ2(w1,w3)+λ2(w2))x + V9(w1,w2,w3)e(µ1(w2,w3)+λ1(w1))x +

V12(w1,w2,w3)e(µ2(w2,w3)+λ2(w1))x + +

V15(w1,w2,w3)e

+

(λ1(w3)+λ2(w1)+λ2(w2))x

V16(w1,w2,w3)e

+

(λ1(w1)+λ1(w2)+λ2(w3))x

V17(w1,w2,w3)e

+

V18(w1,w2,w3)e

+

(λ1(w1)+λ2(w2)+λ2(w3))x

V20(w1,w2,w3)e

V10 )

2Vcb˜ ccU1,1D2,12 3µ2,23λ1,1

V11 )

2Vcb˜ ccU1,1D2,23 3µ1,23λ2,1

V12 )

2Vcb˜ ccU2,1D2,23 3µ2,23λ2,1

V13 ) Vc[2/3b˜ cc(d1,13U1,2 + U1,3d1,12 + U1,1d1,23) + c˜cccU1,1U1,2U1,3] λ1,1λ1,2 + λ1,3λ1,2 + λ1,1λ1,3

Vc[2/3b˜ cc(d3,12U1,3 + U1,2d3,13 + U2,1d1,23) + c˜cccU2,1U1,2U1,3] λ2,1λ1,2 + λ1,3λ1,2 + λ2,1λ1,3 V16 ) Vc[2/3b˜ cc(d4,12U1,3 + U2,2d3,13 + U2,1d3,23 + c˜cccU2,1U2,2U1,3)] λ2,1λ2,2 + λ1,3λ2,2 + λ2,1λ1,3 V17 )

V19(w1,w2,w3)e(λ1(w2)+λ2(w1)+λ2(w3))x + (λ2(w1)+λ2(w2)+λ2(w3))x

2Vcb˜ ccU1,1D1,23 3µ1,23λ1,1

V15 )

V13(w1,w2,w3)e(λ1(w1)+λ1(w2)+λ1(w3))x + (λ1(w2)+λ2(w1)+λ1(w3))x

V9 )

Vc[2/3b˜ cc(d2,12U1,3 + U2,2d1,13 + U1,1d3,13) + c˜cccU1,1U2,2U1,3] λ1,1 ‚λ2,2 + λ1,3λ2,2 + λ1,1λ1,3

V11(w1,w2,w3)e(µ1(w2,w3)+λ2(w1))x +

V14(w1,w2,w3)e

2Vcb˜ ccU2,2D1,12 3µ2,13λ2,2

V14 )

V10(w1,w2,w3)e(µ2(w2,w3)+λ1(w1))x +

(λ1(w1)+λ2(w2)+λ1(w3))x

V8 )

(A-22)

Vc[2/3b˜ cc(d1,13U2,3 + U1,2d2,12 + U1,1d2,23) + c˜cccU1,1U1,2U2,3] λ1,1λ1,2 + λ2,3λ1,2 + λ1,1λ2,3 V18 )

In this equation,

χ1,2(w1,w2,w3) ) NTP ( xNTP2 + 2NTPj(1 + a˜cf)(w1 + w2 + w3) (A-23) are the characteristic values of eq A-19, V1(w1,w2,w3) through V20(w1,w2,w3) are the integration constants corresponding to its particular solution:

V1 )

2Vcb˜ ccU1,3D1,12 3µ1,12λ1,3

V2 )

2Vcb˜ ccU1,3D2,12 3µ2,12λ1,3

Vc[2/3b˜ cc(d2,12U2,3 + U2,2d2,13 + U1,1d4,23) + c˜cccU1,1U2,2U2,3] λ1,1λ2,2 + λ2,3λ2,2 + λ1,1λ2,3 V19 ) Vc[2/3b˜ cc(d3,12U2,3 + U1,2d4,23 + U2,1d2,23) + c˜cccU2,1U1,2U2,3] λ1,2λ2,1 + λ2,1λ2,3 + λ1,2λ2,3 V20 ) Vc[2/3b˜ cc(d4,12U2,3 + U2,2d4,23 + U2,1d4,23) + c˜cccU2,1U2,2U2,3] λ2,2λ2,1 + λ2,1λ2,3 + λ2,2λ2,3 and V1(w1,w2,w3) and V2(w1,w2,w3) are the integration constants corresponding to its complementary solution, which are obtained from the following set of algebraic equations:

V4 )

2Vcb˜ ccU2,3D2,12 3µ2,12λ2,3

V1 + V2 + V1 + V2 + V3 + V4 + V5 + ... + V20 ) 1 [χ V + χ2V2 + V1(µ1,12 + λ1,3) + V2(µ2,12 + λ1,3) + 2NTP 1 1 V3(λ2,3 + µ1,12) + V4(λ2,3 + µ2,12) + V5(µ1,13 + λ1,2) + ... + V19(λ1,2 + λ2,1 + λ2,3) + V20(λ2,1 + λ2,2 + λ2,3)] (A-24)

V5 )

2Vcb˜ ccU1,2D1,13 3µ1,13λ1,2

V1χ1eχ1 + V2χ2eχ2 + V1(µ1,12 + λ1,3)e(µ1,12+λ1,3) +

2Vcb˜ ccU2,3D1,12 V3 ) 3µ1,12λ2,3

2Vcb˜ ccU1,2D2,13 V6 ) 3µ2,13λ1,2 2Vcb˜ ccU2,2D1,13 V7 ) 3µ1,13λ2,2

V2(µ2,12 + λ1,3)e(µ2,12+λ1,3) + V3(µ1,12 + λ2,3)e(µ1,12+λ2,3) + V4(µ2,12 + λ2,3)e(µ2,12+λ2,3) + V5(µ1,12 + λ1,3)e(µ1,12+λ2,3) + ... + V19(λ1,2 + λ2,1 + λ2,3)e(λ1,2+λ2,1+λ2,3) + V20(λ2,1 + λ2,2 + λ2,3)e(λ2,1+λ2,2+λ2,3) ) 0 (A-25)


The following shortened notations have been introduced here: Vi ) Vi(w1,w2,w3), Vi ) Vi(w1,w2,w3), χi ) χi(w1,w2,w3), Vc ) NTPfj(w1 + w2 + w3), Di,jk ) Di(wj,wk), di,jk ) di(wj,wk), and µi,jk ) µi(wj,wk). For x ) 1, the function F3,ccc(w1,w2,w3) (eq A-22) reduces to the function G3,ccc(w1,w2,w3), which is given by eq 40 in Section 3.2.1. A2. The Gθ Functions. A2.1. The First-Order Gθ Function. This function is obtained by collecting the terms containing B1eju1t in the equations obtained in Step 4, and equating them to zero. The following ODE and its boundary conditions are obtained:

d2F1,θ(x,u1) dx2

dF1,θ(x,u1) - 2NTP dx 2NTPju1F1,θ(x,u1)(1 + a˜cf) ) 2NTPju1a˜θ (A-26)

x ) 0: F1,θ(0,u1) ) x ) 1:

1 dF1,θ(x,u1) | | 2NTP dx |x)0

F2,θθ(x,u1,u2) ) S1(u1,u2)eµ1(u1,u2)x + S2(u1,u2)eµ2(u1,u2)x + s1(u1,u2)e(λ1(u1)+λ1(u2))x + s2(u1,u2)e(λ1(u1)+λ2(u2))x + s3(u1,u2)e(λ2(u1)+λ1(u2))x + s4(u1,u2)e(λ2(u2)+λ2(u1))x + s5(u1,u2)eλ1(u1)x + s6(u1,u2)eλ2(u2)x + s7(u1,u2)eλ1(u2)x + s8(u1,u2)eλ2(u1)x + s9 (A-36) The integration constants s1(u1,u2) through s9 correspond to the particular solution and can be calculated using the following expressions:

s1 ) s2 )

(A-27) s3 )

dF1,θ(x,u1) | | )0 dx |x)1

(A-28) s4 )

NTP j(u1 + u2) f b˜ ccH1,1H1,2 λˆ 1,1λˆ 1,2 NTP j(u1 + u2) f b˜ ccH1,1H2,2 λˆ 2,2λˆ 1,1 NTP j(u1 + u2) f b˜ ccH1,2H2,1 λˆ 2,1λˆ 1,2 NTP j(u1 + u2) f b˜ ccH2,1H2,2

The solution is

F1,θ(x,u1) ) H1(u1)eλ1(u1)x + H2(u1)eλ2(u1)x + h

(A-29)

s5 ) -

where

h)-

s6 ) -

f a˜θ 1 + a˜c f

(A-30) s7 ) -

and the integration constants H1(u1) and H2(u1) are obtained from the following set of algebaic equations:

H1(u1) + H2(u2) + h )

(H1(u1)λ1(u1) + H2(u1)λ2(u1))/2 NTP (A-31)

H1(u1)λ1(u1)eλ1(u1) + H2(u1)λ2(u1)eλ2(u1) ) 0 (A-32) For x ) 1, eq A-29 reduces to eq 41 in Section 3.2.2, defining G1,θ(u1). A2.2. The Second-Order Gθ Function. By collecting the terms that contain B1B2ej(u1+u2)t in the equations obtained in Step 4 and by equating them to zero, the following is obtained:

d2F2,θθ(x,u1,u2)

dF2,θθ(x,u1,u2) - 2NTP dx dx 2NTP j(u1 + u2)F2,θθ(x,u1,u2)(1 + a˜c f) ) 2

[

]

(F1,θ(x,u1) + F1,θ(x,u2)) (A-33) 2

1 dF2,θθ(x,u1,u2) | (A-34) x ) 0: F2,θθ(0,u1,u2) ) | 2NTP dx |x)0 x ) 1: The solution is

dF2,θθ(x,u1,u2) | | )0 dx |x)1

(u1 + u2) f H1,1[(b˜ cθ/2) - hb˜ cc] u2(1 + a˜c f ) (u1 + u2) f H2,2[(b˜ cθ/2) - hb˜ cc] u1(1 + a˜c f ) (u1 + u2) f H1,2[(b˜ cθ/2) - hb˜ cc] u1(1 + a˜c f ) (u1 + u2) f H2,1[(b˜ cθ/2) - hb˜ cc]

s9 ) -

u2(1 + a˜c f) f (b˜ θθ + b˜ cch2 - b˜ cθh) (1 + a˜c f )

whereas the integration constants S1(u1,u2) and S2(u1,u2), which correspond to the complementary solution, are obtained from the following set of algebraic equations:

S 1 + S 2 + s1 + s2 + s3 + s4 + s 5 + s6 + s7 + s8 + s 9 ) {[µˆ 1,12S1 + µˆ 2,12S2 + s1(λˆ 1,1 + λˆ 1,2) + s2(λˆ 1,1 + λˆ 2,2) + s3(λˆ 2,1 + λˆ 1,2) + s4(λˆ 2,1 + λˆ 2,2) + s5λˆ 1,1 + s6λˆ 2,2 + s7λˆ 1,2 + s8λˆ 2,1]/2}/NTP (A-37) S1µˆ 1,12eµˆ 1,12 + S2µˆ 2,12eµˆ 2,12 + s1(λˆ 1,1 + λˆ 2,2)e(λˆ 1,1+λˆ 1,2) +

2NTP j(u1 + u2)f b˜ θθ + b˜ ccF1,θ(x,u1)F1,θ(x,u2) + b˜ cθ

s8 ) -

λˆ 2,1λˆ 2,2

(A-35)

s2(λˆ 1,1 + λˆ 2,2)e(λˆ 1,1+λˆ 2,2) + s3(λˆ 2,1 + λˆ 1,2)e(λˆ 2,1+λˆ 1,2) + s4(λˆ 2,1 + λˆ 2,2)e(λˆ 2,1+λˆ 2,2) + s5λˆ 1,1eλˆ 1,1 + s6λˆ 2,2eλˆ 2,2 + s7λˆ 1,2eλˆ 1,2 + s8λˆ 2,1 ) 0 (A-38) The following shortened notations have been introduced here: si ) si(u1,u2), Si ) Si(u1,u2), Hi,j ) Hi(uj), λˆ i,j ) λi(uj), and µˆ i,jk ) µi(uj,uk). For x ) 1 in eq A-36, the final expression for G2,θθ(u1,u2), which is given by eq 42 in Section 3.2.2, is obtained. A2.3. The Third-Order Gθ Function. The third-order Gθ function is obtained by collecting the terms of the equations obtained in Step 4 that contain B1B2B3ej(u1+u2+u3)t and by equating

284


them to zero. The following ODE and its boundary conditions are obtained:

The integration constants p1(u1,u2,u3) through p45(u1,u2,u3) are defined in the following way:

d2F3,θθθ(x,u1,u2,u3)

dF3,θθθ(x,u1,u2,u3) - 2NTP dx dx 2NTPj(u1 + u2 + u3)F3,θθθ(x,u1,u2,u3)(1 + a˜cf) ) 2

[

2

p1 )

/3b˜ cc(F2,θθ(x,u1,u2)F1,θ(x,u3) + F2,θθ(x,u1,u3)F1,θ(x,u2) +

F2,θθ(x,u2,u3)F1,θ(x,u1)) + 1/3b˜ cθ(F2,θθ(u1,u2) + F2,θθ(u1,u3) + F2,θθ(u2,u3)) + c˜cccF1,θ(x,u1)F1,θ(x,u2)F1,θ(x,u3) + c˜θθθ + c˜ccθ (F (x,u1)F1,θ(x,u2) + F1,θ(x,u1)F1,θ(x,u3) + 3 1,θ c˜cθθ (F (x,u1) + F1,θ(x,u2) + F1,θ(x,u2)F1,θ(x,u3)) + 3 1,θ

]

F1,θ(x,u3)) 2NTPj(u1 + u2 + u3)f (A-39) x ) 0: F3,θθθ(0,u1,u2,u3) )

1 dF3,θθθ(x,u1,u2,u3) | (A-40) | 2NTP dx |x)0

dF3,θθθ(x,u1,u2,u3) | x ) 1: | )0 dx |x)1

p2 )

p3 )

p4 )

p5 )

(A-41) p6 )

Its solution is F3,θθθ(x,u1,u2,u3) ) P1(u1,u2,u3)eχ1(u1,u2,u3)x +

p7 )

P2(u1,u2,u3)eχ2(u1,u2,u3)x + p1(u1,u2,u3)e(µ1(u1,u2)+λ1(u3))x + p2(u1,u2,u3)e

(µ2(u1,u2)+λ1(u3))x

p4(u1,u2,u3)e

(µ2(u1,u2)+λ2(u3))x

+ p3(u1,u2,u3)e

(µ1(u1,u2)+λ2(u3))x

+

+ p5(u1,u2,u3)e

(µ1(u1,u3)+λ1(u2))x

+

p8 )

p6(u1,u2,u3)e(µ2(u1,u3)+λ1(u2))x + p7(u1,u2,u3)e(µ1(u1,u3)+λ2(u2))x + p8(u1,u2,u3)e(µ2(u1,u3)+λ2(u2))x + p9(u1,u2,u3)e(µ1(u2,u3)+λ1(u1))x + (µ2(u2,u3)+λ1(u1))x

p10(u1,u2,u3)e

+ p11(u1,u2,u3)e

(µ1(u2,u3)+λ2(u1))x

+

p12(u1,u2,u3)e(µ2(u2,u3)+λ2(u1))x + p13(u1,u2,u3)eµ1(u1,u2)x + p14(u1,u2,u3)eµ2(u1,u2)x + p15(u1,u2,u3)eµ1(u1,u3)x +

p9 )

p10 )

p16(u1,u2,u3)eµ2(u1,u3)x + p17(u1,u2,u3)eµ1(u2,u3)x + p18(u1,u2,u3)eµ2(u2,u3)x + p19(u1,u2,u3)e(λ1(u1)+λ1(u2) + λ1(u3))x + p20(u1,u2,u3)e(λ1(u1)+λ2(u2)+λ1(u3))x + (λ1(u2)+λ2(u1)+λ1(u3))x

+

(λ1(u3)+λ2(u1)+λ2(u2))x

+

p21(u1,u2,u3)e p22(u1,u2,u3)e

p11 )

p12 )

p23(u1,u2,u3)e(λ1(u1)+λ1(u2)+λ2(u3))x + p24(u1,u2,u3)e(λ1(u1)+λ2(u2)+λ2(u3))x +

p13 ) -

p25(u1,u2,u3)e(λ1(u2)+λ2(u1)+λ2(u3))x + p26(u1,u2,u3)e(λ2(u1)+λ2(u2)+λ2(u3))x + p27(u1,u2,u3)e(λ1(u1)+λ1(u3))x + p28(u1,u2,u3)e(λ1(u1)+λ1(u2))x + p29(u1,u2,u3)e(λ1(u2)+λ1(u3))x + p30(u1,u2,u3)e(λ2(u2)+λ1(u3))x + p31(u1,u2,u3)e(λ2(u1)+λ1(u3))x + p32(u1,u2,u3)e(λ1(u1)+λ2(u3))x + p33(u1,u2,u3)e(λ2(u2)+λ2(u3))x + p34(u1,u2,u3)e(λ1(u2)+λ2(u3))x + p35(u1,u2,u3)e(λ2(u1)+λ2(u3))x + p36(u1,u2,u3)e(λ1(u1)+λ2(u2))x + p37(u1,u2,u3)e(λ2(u1)+λ1(u2))x + p38(u1,u2,u3)e(λ2(u1)+λ2(u2))x + λ1(u1)x λ1(u3)x

p39(u1,u2,u3)e p41(u1,u2,u3)e p43(u1,u2,u3)e

λ2(u2)x

+ p40(u1,u2,u3)e

λ1(u2)x

+

+ p42(u1,u2,u3)e

λ2(u1)x

+

+ p44(u1,u2,u3)e

λ2(u3)x

+ p45(u1,u2,u3) (A-42)

p14 ) p15 ) p16 ) p17 ) p18 ) -

pc1S1,12H1,3 2µˆ 1,12λˆ 1,3 pc1S2,12H1,3 2µˆ 2,12λˆ 1,3 pc1S1,12H2,3 2µˆ 1,12λˆ 2,3 pc1S2,12H2,3 2µˆ 2,12λˆ 2,3 pc1S1,13H1,2 2µˆ 1,13λˆ 1,2 pc1S2,13H1,2 2µˆ 2,13λˆ 1,2 pc1S1,13H2,2 2µˆ 1,13λˆ 2,2 pc1S2,13H2,2 2µˆ 2,13λˆ 2,2 pc1S1,23H1,1 2µˆ 1,23λˆ 1,1 pc1S2,23H1,1 2µˆ 2,23λˆ 1,1 pc1S1,23H2,1 2µˆ 1,23λˆ 2,1 pc1S2,23H2,1 2µˆ 2,23λˆ 2,1 pc2S1,12

2NTP j(1 + a˜cf)u3 pc2S2,12 2NTP j(1 + a˜cf)u3 pc2S1,13 2NTP j(1 + a˜cf)u2 pc2S2,13 2NTP j(1 + a˜cf)u2 pc2S1,23 2NTP j(1 + a˜cf)u1 pc2S2,23 2NTP j(1 + a˜cf)u1


p19 )

p20 )

p21 )

p22 )

p23 )

p24 )

p25 )

p26 )

p27 )

p28 )

p29 )

p30 )

p31 )

p32 ) p33 )

p34 )

p35 )

p36 )

pc3[2/3b˜ cc(s1,12H1,3 + H1,2s1,13 + H1,1s1,23) + c˜cccH1,1H1,2H1,3] λˆ 1,1λˆ 1,2 + λˆ 1,3λˆ 1,2 + λˆ 1,1λˆ 1,3 pc3[2/3b˜ cc(s2,12H1,3 + H2,2s1,13 + H1,1s3,23) + c˜cccH1,1H2,2H1,3] λˆ 1,1λˆ 2,2 + λˆ 1,3λˆ 2,2 + λˆ 1,1λˆ 1,3 pc3[2/3b˜ cc(s3,12H1,3 + H1,2s3,13 + H2,1s1,23) + c˜cccH2,1H1,2H1,3] λˆ 2,1λˆ 1,2 + λˆ 1,3λˆ 1,2 + λˆ 2,1λˆ 1,3 pc3[2/3b˜ cc(s4,12H1,3 + H2,2s3,13 + H2,1s3,23) + c˜cccH2,1H2,2H1,3] λˆ 2,1λˆ 2,2 + λˆ 1,3λˆ 2,2 + λˆ 2,1λˆ 1,3 pc3[2/3b˜ cc(s1,12H2,3 + H1,2s2,13 + H1,1s2,23) + c˜cccH1,1H1,2H2,3] λˆ 1,1λˆ 1,2 + λˆ 2,3λˆ 1,2 + λˆ 1,1λˆ 2,3 pc3[2/3b˜ cc(s2,12H2,3 + H2,2s2,13 + H1,1s4,23) + c˜cccH1,1H2,2H2,3] λˆ 1,1λˆ 2,2 + λˆ 2,3λˆ 2,2 + λˆ 1,1λˆ 2,3 pc3[2/3b˜ cc(s3,12H2,3 + H2,1s2,23 + H1,2s4,13) + c˜cccH2,1H1,2H2,3] λˆ 1,2λˆ 2,1 + λˆ 2,1λˆ 2,3 + λˆ 1,2λˆ 2,3 pc3[2/3b˜ cc(s4,12H2,3 + H2,2s4,13 + H2,1s4,23) + c˜cccH2,1H2,2H2,3] λˆ 2,2λˆ 2,1 + λˆ 2,1λˆ 2,3 + λˆ 2,2λˆ 2,3

pc3{2/3b˜ cc(s5,12H1,3 - hs1,13 + H1,1s7,23) + 1/3b˜ cθs1,13 + [(c˜ccθ/3) - hc˜ccc]H1,1H1,3} λˆ 1,1λˆ 1,3 - NTP j(1 + a˜cf)u2 pc3{2/3b˜ cc(s5,13H1,2 - hs1,12 + H1,1s5,23) + 1/3b˜ cθs1,12 + [(c˜ccθ/3) - hc˜ ccc]H1,1H1,2} λˆ 1,1λˆ 1,2 - NTP j(1 + a˜cf)u3 pc3{2/3b˜ cc(s7,12H1,3 - hs1,23 + H1,2s7,13) + 1/3b˜ cθs1,23 + [(c˜ccθ/3) - hc˜ccc]H1,2H1,3} λˆ 1,2λˆ 1,3 - NTP j(1 + a˜cf)u1 pc3{2/3b˜ cc(s6,12H1,3 - hs3,23 + H2,2s7,13) + 1/3b˜ cθs3,23 + [(c˜ccθ/3) - hc˜ccc]H2,2H1,3} λˆ 2,2λˆ 1,3 - NTP j(1 + a˜cf)u1 pc3{2/3b˜ cc(s8,12H1,3 - hs3,13 + H2,1s7,23) + 1/3b˜ cθs3,13 + [(c˜ccθ/3) - hc˜ccc]H2,1H1,3} λˆ 2,1λˆ 1,3 - NTP j(1 + a˜cf)u2 pc3{2/3b˜ cc(s5,12H2,3 - hs2,13 + H1,1s6,23) + 1/3b˜ cθs2,13 + [(c˜ccθ/3) - hc˜ccc]H1,1H2,3} λˆ 1,1λˆ 2,3 - NTP j(1 + a˜cf)u2 pc3{2/3b˜ cc(s6,12H2,3 - hs4,23 + H2,2s6,13) + 1/3b˜ cθs4,23 + [(c˜ccθ/3) - hc˜ccc]H2,3H2,2} λˆ 2,2λˆ 2,3 - NTP j(1 + a˜cf)u1 pc3{2/3b˜ cc(s7,12H2,3 - hs2,23 + H1,2s6,15) + 1/3b˜ cθs2,23 + [(c˜ccθ/3) - hc˜ccc]H2,3H1,2} λˆ 1,2λˆ 2,3 - NTP j(1 + a˜cf)u1 pc3{2/3b˜ cc(s8,12H2,3 - hs4,13 + H2,1s6,23) + 1/3b˜ cθs4,13 + [(c˜ccθ/3) - hc˜ccc]H2,1H2,3} λˆ 2,1λˆ 2,3 - NTP j(1 + a˜cf)u2 pc3{2/3b˜ cc(s8,23H1,1 - hs2,12 + H2,2s5,13) + 1/3b˜ cθs2,12 + [(c˜ccθ/3) - hc˜ccc]H1,1H2,2} λˆ 1,1λˆ 2,2 - NTP j(1 + a˜cf)u3

286


p37 )

p38 ) p39 ) p40 ) p41 ) p42 ) p43 ) p44 ) -

pc3{2/3b˜ cc(s8,13H1,2 - hs3,12 + H2,1s5,23) + 1/3b˜ cθs3,12 + [(c˜ccθ/3) - hc˜ccc]H2,1H1,2} λˆ 2,1λˆ 1,2 - NTP j(1 + a˜cf)u3 pc3{2/3b˜ cc(s8,13H2,2 - hs4,12 + H2,1s8,23) + 1/3b˜ cθs4,12 + [(c˜ccθ/3) - hc˜ccc]H2,2H2,1} λˆ 2,1λˆ 2,2 - NTP j(1 + a˜cf)u3

pc3{2/3b˜ cc(H1,1s9,23 - hs5,12 - hs5,13) + 1/3b˜ cθ(s5,13 + s5,12) + H1,1[3c˜ccch2 - c˜ccθh + (c˜ccθ)]/3} NTP j(1 + a˜cf)(u2 + u3) pc3{2/3b˜ cc(H1,2s9,13 - hs7,12 - hs5,23) + 1/3b˜ cθ(s5,23 + s7,12) + H1,2[3c˜ccch2 - c˜ccθh + (c˜ccθ)]/3} NTP j(1 + a˜cf)(u1 + u3) pc3{2/3b˜ cc(H1,3s9,12 - hs7,13 - hs7,23) + 1/3b˜ cθ(s7,13 + s7,23) + H1,3[3c˜ccch2 - c˜ccθh + (c˜cθθ)]/3} NTP j(1 + a˜cf)(u1 + u2) pc3{2/3b˜ cc(H2,1s9,23 - hs8,12 - hs8,13) + 1/3b˜ cθ(s8,12 + s8,13) + H2,1[3c˜ccch2 - c˜ccθh + (c˜cθθ)]/3} NTP j(1 + a˜cf)(u2 + u3) pc3{2/3b˜ cc(H2,2s9,13 - hs6,12 - hs8,23) + 1/3b˜ cθ(s6,12 + s8,23) + H2,2[3c˜ccch2 - c˜ccθh + (c˜cθθ)]/3} NTP j(1 + a˜cf)(u1 + u3) pc3{2/3b˜ cc(H2,3s9,12 - hs6,13 - hs6,23) + 1/3b˜ cθ(s6,13 + s6,23) + H2,3[3c˜ccch2 - c˜ccθh + (c˜cθθ)]/3} NTP j(1 + a˜cf)(u1 + u2) p45 ) -

f{(s9,12 + s9,13 + s9,23)(b˜ cθ - 2b˜ cch)/3 + c˜ccθh2 - c˜ccch3 + c˜θθθ - c˜cθθh} (1 + a˜cf)

The integration constants P1(u1,u2,u3) and P2(u1,u2,u3) are obtained by solving the following set of algebraic equations:

P1 + P2 + p1 + p2 + p3 + p4 + p5 + p6 + ... + p45 ) {[χˆ 1P1 + χˆ 2P2 + p1(µˆ 1,12 + λˆ 1,3) + p2(µˆ 2,12 + λˆ 1,3) + p3(λˆ 2,3 + µˆ 1,12) + p4(λˆ 2,3 + µˆ 2,12) + p5(µˆ 1,13 + λˆ 1,2) + ... + p43λˆ 2,2 + p44λˆ 2,3]/2}/NTP (A-43)

d2F2,cθ(x,w1,u1)

dF2,cθ(x,w1,u1) dx dx 2NTP j(w1 + u1)F2,cθ(x,w1,u1)(1 + a˜c f) ) 2NTP j(w1 + u1) × f [b˜ ccF1,c(x,w1)F1,θ(x,u1) + b˜ cθ/2F1,c(x,w1)] (A-45) 2

- 2NTP

x ) 0: F2,cθ(0,w1,u1) )

P1χˆ 1eχˆ 1 + P2χˆ 2eχˆ 2 + p1(µˆ 1,12 + λˆ 1,3)e(µˆ 1,12+λˆ 1,3) + p2(µˆ 2,12 + λˆ 1,3)e(µˆ 2,12+λˆ 1,3) + p3(µˆ 1,12 + λˆ 2,3)e(µˆ 1,12+λˆ 2,3) + p4(µˆ 2,12 + λˆ 2,3)e

(µˆ 2,12+λˆ 2,3)

+ p5(µˆ 1,13 + λˆ 1,2)e p43λˆ 2,2e

λˆ 2,2

(µˆ 1,13+λˆ 1,2) λˆ 2,3

+ p44λˆ 2,3e

x ) 1:

+ ... +

) 0 (A-44)

The following shortened notations were introduced here: pi ) pi(u1,u2,u3), Pi ) Pi(u1,u2,u3), pc1 ) 4/3NTPfj(u1 + u2 + u3)b˜ cc, pc2 ) 2NTPfj(u1 + u2 + u3)(b˜ cθ - 2b˜ cch)/3, pc3 ) NTPfj(u1 + u2 + u3), Si,jk ) Si(uj,uk), si,jk ) si(uj,uk), and χˆ i ) χi(uj,uk,ul). The function G3,θθθ(u1,u2,u3) (eq 43 in Section 3.2.2) is identical to the function F3,θθθ(x,u1,u2,u3) (given in eq A-42) for x ) 1. A3. The Gcθ Functions. A3.1. The Second-Order Gcθ Function. To derive the second-order Gcθ function, G2,cθ(w1,u1), it is necessary to collect the terms of the equations obtained in Step 4 that contain A1B2ej(w1+u2)t and equate them to zero. The resulting nonhomogeneous ODE and its boundary conditions are

1 dF2,cθ(x,w1,u1) | | (A-46) 2NTP dx |x)0

dF2,cθ(x,w1,u1) | | )0 dx |x)1

(A-47)

The solution is

F2,cθ(x,w1,u1) ) K1(w1,u1)eµ1(w1,u1)x + K2(w1,u1)eµ2(w1,u1)x + k1(w1,u1)e(λ1(w1)+λ1(u1))x + k2(w1,u1)e(λ1(w1)+λ2(u1))x + k3(w1,u1)e(λ2(w1)+λ1(u1))x + k4(w1,u1)e(λ2(w1)+λ2(u1))x + k5(w1,u1)eλ1(w1)x + k6(w1,u1)eλ2(w1)x (A-48) with the following integration constants k1(w1,u1)-k6(w1,u1):

k1 )

k2 )

NTP j(w1 + u1) fb˜ ccU1,1H1,1 λ1,1λˆ 1,1 NTP j(w1 + u1) fb˜ ccU1,1H2,1 λ1,1λˆ 2,1


NTP j(w1 + u1) fb˜ ccU2,1H1,1 λˆ 1,1λ2,1

k3 )

F3,ccθ(x,w1,w2,u1) ) Z1(w1,w2,u1)eχ1(w1,w2,u1)x + Z2(w1,w2,u1)eχ2(w1,w2,u1)x + z1(w1,w2,u1)e(µ1(w1,u1)+λ1(w2))x + z2(w1,w2,u1)e(µ2(w1,u1)+λ1(w2))x + z3(w1,w2,u1)e(µ1(w1,u1)+λ2(w2))x +

NTP j(w1 + u1) fb˜ ccU2,1H2,1 k4 ) λˆ 2,1λ2,1 k5 ) k6 ) -

z4(w1,w2,u1)e(µ2(w1,u1)+λ2(w2))x + z5(w1,w2,u1)e(µ1(w2,u1)+λ1(w1))x + z6(w1,w2,u1)e(µ2(w2,u1)+λ1(w1))x + z7(w1,w2,u1)e(µ1(w2,u1)+λ2(w1))x +

(w1 + u1)fU1,1[(b˜ cθ/2) - hb˜ cc]

z8(w1,w2,u1)e(µ2(w2,u1)+λ2(w1))x + z9(w1,w2,u1)e(µ1(w1,w2)+λ1(u1))x +

u1(1 + a˜c f)

z10(w1,w2,u1)e(µ2(w1,w2)+λ1(u1))x +

(w1 + u1)fU2,1[(b˜ cθ/2) - hb˜ cc]

z11(w1,w2,u1)e(µ1(w1,w2)+λ2(u1))x +

u1(1 + a˜c f)

z12(w1,w2,u1)e(µ2(w1,w2)+λ2(u1))x + z13(w1,w2,u1)eµ1(w1,w2)x +

and K1(w1,u1) and K2(w1,u1) obtained as solutions of the following set of equations:

{[µˇ 1,11K1 + µˇ 2,11K2 + k1(λ1,1 + λˆ 1,1) + k2(λ1,1 + λˆ 2,1) + k3(λ2,1 + λˆ 1,1) + k4(λ2,1 + λˆ 2,2) + k5λ1,1 + k6λ2,1]/2}/NTP (A-49) K1µˇ 1,11eµˇ 1,11 + K1µˇ 1,11eµˇ 1,11 + k1(λ1,1 + λˆ 1,1)e(λ1,1+λˆ 1,1) +

z18(w1,w2,u1)e(λ1(w2)+λ2(w1)+λ2(u1))x + z19(w1,w2,u1)e(λ1(w1)+λ2(w2)+λ1(u1))x + z20(w1,w2,u1)e(λ1(w1)+λ2(w2)+λ2(u1))x + z21(w1,w2,u1)e(λ2(w1)+λ2(w2)+λ1(u1))x + z22(w1,w2,u1)e(λ2(w1)+λ2(w2)+λ2(u1))x +

k2(λ1,1 + λˆ 2,1)e(λ1,1+λˆ 2,1) + k3(λ2,1 + λˆ 1,1)e(λ2,1+λˆ 1,1) + k4(λ2,1 + λˆ 2,1)e(λ2,1+λˆ 2,1) + k5λ1,1eλ1,1 + k6λ2,1eλ2,1 ) 0 (A-50) The shortened notations introduced here are as follows: ki ) ki(w1,u1), Ki ) Ki(w1,u1), and µˇ i,jk ) µi(wi,uk). For x ) 1, eq A-48 reduces to eq 44 in Section 3.2.3; i.e., function F2,cθ(x,w1,u1) becomes identical to G2,cθ(w1,u1). A3.2. The Third-Order Gcθ Functions. Two third-order Gcθ functions can be defined: G3,ccθ(w1,w2,u1) and G3,cθθ(w1,u1,u2). A3.2.1. Derivation of G3,ccθ(w1,w2,u1). To derive the function G3,ccθ(w1,w2,u1), it is necessary to collect the terms in the equations obtained in Step 4 that contain A1A2B1ej(w1+w2+u1)t and equate them to zero. The resulting ODE and its boundary conditions are

2

z16(w1,w2,u1)e(λ1(w1)+λ1(w2)+λ2(u1))x + z17(w1,w2,u1)e(λ2(w1)+λ1(w2)+λ1(u1))x +

K 1 + K2 + k1 + k2 + k3 + k4 + k5 + k6 )

d2F3,ccθ(x,w1,w2,u1)

z14(w1,w2,u1)eµ2(w1,w2)x + z15(w1,w2,u1)e(λ1(w1)+λ1(w2)+λ1(u1))x +

dF3,ccθ(x,w1,w2,u1) - 2NTP dx

dx 2NTP j(w1 + w2 + u1)F3,ccθ(x,w1,w2,u1)(1 + a˜c f) ) 2NTP j(w1 + w2 + u1) f{2/3b˜ cc[F2,cθ(x,w1,u1)F1,c(x,w2) +

z23(w1,w2,u1)e(λ1(w1)+λ1(w2))x + z24(w1,w2,u1)e(λ1(w2)+λ2(w1))x + z25(w1,w2,u1)e(λ1(w1)+λ2(w2)) + z26(w1,w2,u1)e(λ2(w1)+λ2(w2))x (A-54) with the integration constants z1(w1,w2,u1) to z26(w1,w2,u1) defined in the following way:

z1 )

zc1K1,11U1,2 2µˇ 1,11λ1,2

z2 )

zc1K2,11U1,2 2µˇ 2,11λ1,2

z3 )

zc1K1,11U2,2 2µˇ 1,11λ2,2

z4 )

zc1K2,11U2,2 2µ j 2,11λ2,2

z5 )

zc1K1,21U1,1 2µˇ 1,21λ1,1

z6 )

zc1K2,21U1,1 2µˇ 2,21λ1,1

z7 )

zc1K1,21U2,1 2µˇ 1,21λ2,1

z8 )

zc1K2,21U2,1 2µˇ 2,21λ2,1

F2,cθ(x,w2,u1)F1,c(x,w1) + F2,cc(x,w1,w2)F1,θ(x,u1)] + /3b˜ cθF2,cc(w1,w2) + c˜cccF1,c(x,w1)F1,c(x,w2)F1,θ(x,u1) +

1

1

z9 )

/3c˜cctF1,c(x,w1)F1,c(x,w2)} (A-51)

x ) 0: F3,cct(0,w1,w2,u) ) x ) 1:

The solution is

1 ∂F3,cct(x,w1,w2,u) | (A-52) | 2N ∂x |x)0

∂F3,ccθ(x,w1,w2,u1) | | )0 ∂x |x)1

z11 )

zc1D1,12H1,1 2µ1,12λˆ 1,1

zc1D1,12H2,1 2µ1,12λˆ 2,1

z10 )

z12 )

zc1D2,12H1,1 2µ2,12λˆ 1,1

zc1D2,12H2,1 2µ2,12λˆ 2,1

(A-53) z13 ) -

zc2D1,12 (1 + a˜c f)u1

z14 ) -

zc2D2,12 (1 + a˜c f )u1

288


z15 )

z16 )

z17 )

z18 )

z19 )

z20 )

z21 )

z22 )

z23 ) z24 ) z25 ) z26 )

zc3[2/3b˜ cc(k1,11U1,2 + U1,1k1,12 + H1,1d1,12) + c˜cccU1,1U1,2H1,1] λ1,1λ1,2 + λ1,1λˆ 1,1 + λ1,2λˆ 1,1 zc3[2/3b˜ cc(k2,11U1,2 + U1,1k2,12 + H2,1d1,12) + c˜cccU1,1U1,2H2,1] λ1,1λ1,2 + λ1,1λˆ 2,1 + λ1,2λˆ 2,1 zc3[2/3b˜ cc(k3,11U1,2 + U2,1k1,12 + H1,1d3,12) + c˜cccU2,1U1,2H1,1] λ2,1λ1,2 + λ2,1λˆ 1,1 + λ1,2λˆ 1,1 zc3[2/3b˜ cc(k4,11U1,2 + U2,1k2,12 + H2,1d3,12) + c˜cccU2,1U1,2H2,1] λ2,1λ1,2 + λ2,1λˆ 2,1 + λ1,2λˆ 2,1 zc3[2/3b˜ cc(k1,11U2,2 + U1,1k3,12 + H1,1d2,12) + c˜cccU1,1U2,2H1,1] λ1,1λ2,2 + λ1,1λˆ 1,1 + λ2,2λˆ 1,1 zc3[2/3b˜ cc(k2,11U2,2 + U1,1k4,12 + H1,1d2,12) + c˜cccU1,1U2,2H2,1] λ1,1λ2,2 + λ1,1λˆ 2,1 + λ2,2λˆ 2,1 zc3[2/3b˜ cc(k3,11U2,2 + U2,1k3,12 + H1,1d4,12) + c˜cccU2,1U2,2H1,1] λ2,1λ2,2 + λ2,1λˆ 1,1 + λ2,2λˆ 1,1 zc3[2/3b˜ cc(k4,11U2,2 + U2,1k4,12 + H2,1d4,12) + c˜cccU2,2U2,2H2,1] λ2,1λ2,2 + λ2,1λˆ 2,1 + λ2,2λˆ 2,1

zc3{2/3b˜ cc(k5,11U1,2 - hd1,12 + U1,1k5,23) + 1/3b˜ cθd1,12 + [(c˜ccθ/3) - hc˜ccc]U1,1U1,2} λ1,1λ1,2 - NTP j(1 + a˜c f)u1 zc3{2/3b˜ cc(k6,11U1,2 - hd3,12 + U2,1k5,12) + 1/3b˜ cθd3,12 + [(c˜ccθ/3) - hc˜ccc]U2,1U1,2} λ2,1λ1,2 - NTP j(1 + a˜c f)u1 zc3{2/3b˜ cc(k5,11U2,2 - hd2,12 + U1,1k6,12) + 1/3b˜ cθd2,12 + [(c˜ccθ/3) - hc˜ccc]U1,1U2,2} λ1,1λ2,2 - NTP j(1 + a˜c f)u1 zc3{2/3b˜ cc(k6,11U2,2 - hd4,12 + U2,1k6,12) + 1/3b˜ cθd4,12 + [(c˜ccθ/3) - hc˜ccc]U2,1U2,2} λ2,1λ2,2 - NTP j(1 + a˜c f)u1

and Z1(w1,w2,u1) and Z2(w1,w2,u1) obtained as solutions of the following set of equations:

Z1 + Z2 + z1 + z2 + z3 + z4 + z5 + z6 + ... + z26 ) {[χˇ 1Z1 + χˇ 2Z2 + z1(µˇ 1,11 + λ1,2) + z2(µˇ 2,11 + λ1,2) + z3(λ2,2 + µˇ 1,11) + z4(λ2,2 + µˇ 2,11) + z5(µˇ 1,21 + λ1,1) + ... + z25(λ1,1 + λ2,2) + z26(λ2,1 + λ2,2)]/2}/NTP (A-55) χˇ 1Z1eχˇ 1 + χˇ 2Z2eχˇ 2 + z1(µˇ 1,11 + λ1,2)e(µˇ 1,11+λ1,2) + z2(µˇ 2,11 + λ1,2)e(µˇ 2,11+λ2,2) + z3(λ2,2 + µˇ 1,11)e(µˇ 1,11+λ2,2) + z4(λ2,2 + µˇ 2,11)e(µˇ 2,11+λ2,2) + z5(µˇ 1,21 + λ1,1)e(µˇ 1,21+λ1,1) + ... + z25(λ1,1 + λ2,2)e(λ1,1+λ2,2) + z26(λ2,1 + λ2,2)e(λ2,2+λ2,1) ) 0 (A-56) The new shortened notations introduced here are as follows: zi ) zi(w1,w2,w3), Zi ) Zi(w1,w2,w3), χˇ i ) χi(w1,w2,u1), zc1 ) 4/ N fj(w + w + u )b ˜ cθ - 2b˜ cch)/3, 3 TP 1 2 1 ˜ cc, zc2 ) f(w1 + w2 + u1)(b zc3 ) NTPfj(w1 + w2 + u1), Ki,jk ) Ki(wj,uk), ki,jk ) ki(wj,uk), and µˇ i,jk ) µi(wj,uk). For x ) 1, eq 45 in Section 3.2.3, which defines G3,ccθ(w1,w2,u1), is obtained from eq A-54.

A3.3.2. Derivation of G3,cθθ(w1,u1,u2). To derive the function G3,cθθ(w1,u1,u2), the terms in the equations obtained in Step 4 that contain A1B1B2ej(w1+u1+u2)t are collected and equated to zero. The resulting ODE and its boundary conditions are

d2F3,cθθ(x,w1,u1,u2)

dF3,cθθ(x,w1,u1,u2) dx dx 2NTP j(w1 + u1 + u2)F3,cθθ(x,w1,u1,u2)(1 + a˜c f ) ) 2

- 2NTP

{2/3b˜ cc[F2,cθ(x,w1,u1)F1,θ(x,u2) + F2,cθ(x,w1,u2)F1,θ(x,u1) + F2,θθ(x,u1,u2)F1,c(x,w1)] + 1/3b˜ cθ[F2,cθ(x,w1,u1) + F2,cθ(x,w1,u2)] + c˜ccc[F1,c(x,w1)F1,θ(x,u1)F1,θ(x,u2)] + 1

/3c˜ccθ[F1,c(x,w1)F1,θ(x,u1) + F1,c(x,w1)F1,θ(x,u2)] + c˜cθθ/3F1,c(x,w1)}2NTP j(w1 + u1 + u2) (A-57) 1 dF3,cθθ(x,w1,u1,u2) | | 2NTP dx |x)0 (A-58) dF3,cθθ(x,w1,u1,u2) | | )0 (A-59) x ) 1: dx |x)1

x ) 0: F3,cθθ(0,w1,u1,u2) )


The solution is

r9 ) F3,cθθ(x,w1,u1,u2) ) R1(w1,u1,u2)eχ1(w1,u1,u2)x + R2(w1,u1,u2)e

χ2(w1,u1,u2)x

+ r1(w1,u1,u2)e

r2(w1,u1,u2)e

(µ2(u1,u2)+λ1(w1))x

r4(w1,u1,u2)e

(µ2(u1,u2)+λ2(w1))x

r6(w1,u1,u2)e

(µ2(w1,u2)+λ1(u1))x

(µ1(u1,u2)+λ1(w1))x

+

+ r3(w1,u1,u2)e

(µ1(u1,u2)+λ2(w1))x

+

+ r5(w1,u1,u2)e

(µ1(w1,u2)+λ1(u1))x

+

+ r7(w1,u1,u2)e

(µ1(w1,u2)+λ2(u1))x

+

r10 )

r11 )

r8(w1,u1,u2)e(µ2(w1,u2)+λ2(u1))x + r9(w1,u1,u2)e(µ1(w1,u1)+λ1(u2))x + r10(w1,u1,u2)e(µ2(w1,u1)+λ1(u2))x + r11(w1,u1,u2)e(µ1(w1,u1)+λ2(u2))x +

r12 )

r12(w1,u1,u2)e(µ2(w1,u1)+λ2(u2))x + r13(w1,u1,u2)eµ1(w1,u2)x + r14(w1,u1,u2)eµ2(w1,u2)x + r15(w1,u1,u2)eµ1(w1,u1)x + r16(w1,u1,u2)e

µ2(w1,u1)x

+ r17(w1,u1,u2)e

r18(w1,u1,u2)e

(λ1(w1)+λ1(u1)+λ1(u2))x

(λ1(w1)+λ1(u1)+λ2(u2))x

+

+

r19(w1,u1,u2)e

(λ2(w1)+λ1(u1)+λ1(u2))x

+

r20(w1,u1,u2)e

(λ1(w1)+λ2(u1)+λ2(u2))x

+

(λ2(w1)+λ1(u1)+λ1(u2))x

+

r22(w1,u1,u2)e

(λ2(w1)+λ2(u2)+λ1(u1))x

+

r23(w1,u1,u2)e

(λ2(w1)+λ2(u1)+λ1(u2))x

+

r24(w1,u1,u2)e

(λ2(w1)+λ2(u1)+λ2(u2))x

+

r21(w,u1,u2)e

r14 ) r15 ) r16 ) -

+ r26(w1,u1,u2)e

(λ1(w1)+λ2(u2))x

+

r27(w1,u1,u2)e

(λ1(w1)+λ1(u2))x

+ r28(w1,u1,u2)e

(λ2(w1)+λ2(u1))x

+

r29(w1,u1,u2)e

(λ2(w1)+λ1(u1))x

+ r30(w1,u1,u2)e

(λ2(w1)+λ2(u2))x

+

r25(w1,u1,u2)e

(λ1(w1)+λ1(u1))x

r31(w1,u1,u2)e(λ2(w1)+λ1(u2))x + r32(w1,u1,u2)e(λ2(w1)+λ2(u1))x + r33(w1,u1,u2)eλ1(w1)x + r34(w1,u1,u2)eλ2(w1)x (A-60) with the following definitions of the integration constants r1(w1,u1,u2) through r34(w1,u1,u2):

r1 )

rc1S1,12U1,1 2µˆ 1,12λ1,1

r13 ) -

rc1K1,11H1,2 2µˇ 1,11λˆ 1,2 rc1K2,11H1,2 2µˇ 2,11λˆ 1,2 rc1K1,11H2,2 2µˇ 1,11λˆ 2,2 rc1K2,11H2,2 2µˇ 2,11λˆ 2,2 rc2K1,12

2jNTP(1 + a˜c f)u1 rc2K2,12 2jNTP(1 + a˜c f)u1 rc2K1,11 2jNTP(1 + a˜c f)u2 rc2K2,11 2jNTP(1 + a˜c f)u2

r17 ) rc3[2/3b˜ cc(s1,12U1,1 + H1,1k1,12 + H1,2k1,11) + c˜cccH1,1H1,2U1,1] λˆ 1,1λˆ 1,2 + λˆ 1,1λ1,1 + λˆ 1,2λ1,1 r18 ) rc3[2/3b˜ cc(s2,12U1,1 + H1,1k2,12 + H2,2k1,11) + c˜cccH1,1H2,2U1,1] λˆ 1,1λˆ 2,2 + λˆ 1,1λ1,1 + λˆ 2,2λ1,1 r19 ) rc3[2/3b˜ cc(s3,12U1,1 + H2,1k1,12 + H1,2k2,11) + c˜cccH2,1H1,2U1,1] λˆ 2,1λˆ 1,2 + λˆ 2,1λ1,1 + λˆ 1,2λ1,1 r20 )

rc1S2,12U1,1 r2 ) 2µˆ 2,12λ1,1

rc3[2/3b˜ cc(s4,12U1,1 + H2,1k2,12 + H2,2k2,11) + c˜cccH2,1H2,2U1,1]

rc1S1,12U2,1 r3 ) 2µˆ 1,12λ2,1

r21 )

rc1S2,12U2,1 r4 ) 2µˆ 2,12λ2,1 r5 )

r6 )

r7 )

r8 )

rc1K1,12H1,1 2µˇ 1,12λˆ 1,1 rc1K2,12H1,1 2µˇ 2,12λˆ 1,1

λˆ 2,1λˆ 2,2 + λˆ 2,1λ1,1 + λˆ 2,2λ1,1 rc3[2/3b˜ cc(s1,12U2,1 + H1,1k3,12 + H1,2k3,11) + c˜cccH1,1H1,2U2,1] λˆ 1,1λˆ 1,2 + λˆ 1,1λ2,1 + λˆ 1,2λ2,1 r22 ) rc3[2/3b˜ cc(s2,12U2,1 + H1,1k4,12 + H2,2k3,11) + c˜cccH1,1H2,2U2,1] λˆ 1,1λˆ 2,2 + λˆ 1,1λ2,1 + λˆ 2,2λ2,1 r23 ) rc3[2/3b˜ cc(s3,12U2,1 + H2,1k3,12 + H1,2k4,11) + c˜cccH2,1H1,2U2,1]

rc1K1,12H2,1 2µˇ 1,12λˆ 2,1

λˆ 2,1λˆ 1,2 + λˆ 2,1λ2,1 + λˆ 1,2λ2,1 r24 )

rc1K2,12H2,1

rc3[2/3b˜ cc(s4,12U2,1 + H2,1k4,12 + H2,2k4,11) + c˜cccH2,1H2,2U2,1]

2µˇ 2,12λˆ 2,1

λˆ 2,1λˆ 2,2 + λˆ 2,1λ2,1 + λˆ 2,2λ2,1

290


r25 )

r26 )

r27 )

r28 )

r29 )

r30 )

r31 )

r32 )

r33 ) -

r34 ) -

rc3{2/3b˜ cc(s5,12U1,1 - hk1,11 + H1,1k5,12) + b˜ cθk1,11/3 + [1/3(c˜ccθ) - hc˜ccc]H1,1U1,1} λˆ 1,1λ1,1 - NTP j(1 + a˜c f)u2 rc3{2/3b˜ cc(s5,12U1,1 - hk1,11 + H1,1k5,12) + b˜ cθk1,11/3 + [1/3(c˜ccθ) - hc˜ccc]H1,1U1,1} λˆ 2,2λ1,1 - NTP j(1 + a˜c f)u1 rc3{2/3b˜ cc(s7,12U1,1 - hk2,12 + H1,2k5,11) + b˜ cθk1,12/3 + [1/3(c˜ccθ) - hc˜ccc]H1,2U1,1} λˆ 1,2λ1,1 - NTP j(1 + a˜c f)u1 rc3{2/3b˜ cc(s8,12U1,1 - hk2,11 + H2,1k5,12) + b˜ cθk2,12/3 + [1/3(c˜ccθ) - hc˜ccc]H2,1U1,1} λˆ 2,1λ1,1 - NTP j(1 + a˜c f)u2 rc3{2/3b˜ cc(s5,12U2,1 - hk3,11 + H1,1k6,12) + b˜ cθk3,11/3 + [1/3(c˜ccθ) - hc˜ccc]H1,1U2,1} λˆ 1,1λ2,1 - NTP j(1 + a˜c f)u2 rc3{2/3b˜ cc(s6,12U2,1 - hk4,12 + H2,2k6,11) + b˜ cθk4,12/3 + [1/3(c˜ccθ) - hc˜ccc]H2,2U21} λˆ 2,2λ2,1 - NTP j(1 + a˜c f)u1 rc3{2/3b˜ cc(s7,12U2,1 - hk3,12 + H1,2k6,11) + b˜ cθk3,12/3 + [1/3(c˜ccθ) - hc˜ccc]H1,2U2,1} λˆ 1,2λ2,1 - NTP j(1 + a˜c f)u1 rc3{2/3b˜ cc(s8,12U2,1 - hk4,11 + H2,1k6,12) + b˜ cθk4,11/3 + [1/3(c˜ccθ) - hc˜ccc]H2,1U2,1} λˆ 2,1λ2,1 - NTP j(1 + a˜c f)u2

rc3[2/3b˜ cc(U1,1s9 - hk5,12 - hk5,11) + b˜ cθ(k5,12 + k5,11)/3 + U2,1(3c˜ccch2 - 2c˜ccθh + c˜cθθ)/3] NTP j(1 + a˜c f)(u1 + u2) rc3[2/3b˜ cc(U2,1s9 - hk6,12 - hk6,11) + b˜ cθ(k6,12 + k6,11)/3 + U2,1(3c˜ccch2 - 2c˜ccθh + c˜cθθ)/3] NTP j(1 + a˜c f)(u1 + u2)

and the integration constants R1(w1,u2,u3) and R2(w1,u2,u3) are obtained by solving the following set of algebraic equations:

R1 + R2 + r1 + r2 + r3 + r4 + r5 + r6 + ... + r34 ) {[χj1R1 + χj2R2 + r1(µˆ 1,12 + λ1,1) + r2(µˆ 2,12 + λ1,1) + r3(λ2,1 + µˆ 1,12) + r4(λ2,1 + µˆ 2,12) + r5(µˇ 1,12 + λ1,1) + ... + r33λ1,1 + r34λ2,1]/2}/NTP (A-61) R1χj1eχj1 + R2χj2eχj2 + r1(µˆ 1,12 + λ1,1)e(µˆ 1,12+λ1,1) + r2(µˆ 2,12 + λ1,1)e(µˆ 2,12+λ1,1) + r3(µˆ 1,12 + λ2,1)e(µˆ 1,12+λ2,1) + r4(µˆ 2,12 + λ2,1)e(µˆ 2,12+λ2,1) + r5(µˇ 1,12 + λˆ 1,1)e(µˇ 1,12+λˆ 1,1) + ... + r33λ1,1eλ1,1 + t34λ2,1eλ2,1 ) 0 (A-62) The following shortened notations are introduced here: ri ) ri(w1,u1,u2), Ri ) Ri(w1,u1,u2), χji ) χi(w1,u1,u2), rc1 ) 4/3NTPfj(w1 + u1 + u2)b˜ cc, rc2 ) 2NTPfj(w1 + u1 + u2)(b˜ cθ - 2b˜ cch)/3, and rc3 ) NTPfj(w1 + u1 + u2). For x ) 1, eq A-60 is reduced to eq 46 in Section 3.2.3, which defines G3,cθθ(w1,u1,u2). Notations A ) amplitude, general and of the nondimensional inlet concentration change

a˜ ) first-order coefficient of the Taylor series expansion of the adsorption isotherm B ) amplitude, general and of the nondimensional temperature change b0 ) Langmuir isotherm coefficient (dm3/mol) b˜ ) second-order coefficient of the Taylor series expansion of the adsorption isotherm C ) concentration in the fluid phase (mol/dm3) c ) nondimensional concentration in the fluid phase c˜ ) third-order coefficient of the Taylor series expansion of the adsorption isotherm Deff ) effective dispersion coefficient (cm2/s) f ) porosity factor Gn ) nth order frequency response function K ) number of input harmonics k ) temperature coefficient of the Langmuir isotherm (1/K) L ) column length (cm) N ) number of frequency response functions NTP ) number of theoretical plates Q ) concentration in the solid phase (mol/dm3) Q0 ) concentration in the solid phase at maximal coverage (Langmuir isotherm) (mol/dm3) q ) nondimensional concentration in the solid phase T ) temperature (K) t ) nondimensional time u ) nondimensional frequency


V ) interstitial fluid velocity (cm/s) w ) nondimensional frequency x ) input (general), nondimensional axial coordinate y ) output (general) z ) input (general); axial coordinate (cm) 〈 〉 ) average value Greek Symbols δ ) relative error (%) ) porosity Φ ) adsorption isotherm relation θ ) nondimensional temperature τ ) time (s) ω ) nondimensional frequency Subscripts A ) adsorption app ) approximate c ) with regard to concentration D ) desorption i ) inlet o ) outlet s ) stationary state x ) with regard to input x z ) with regard to input z θ ) with regard to temperature AbbreViations FR ) frequency response FRF ) frequency response function HFRF ) higher-order frequency response function Literature Cited (1) Skarstrom, C. W. Method and apparatus for fractionating gaseous mixtures by adsorption. U.S. Patent No. 2,944,627, 1960. (2) Wilhelm, R. H.; Rice, A. W.; Bendelius, A. R. Parametric Pumping: A Dynamic Principle for Separation of Fluid Mixtures. Ind. Eng. Chem. Fundam. 1966, 5, 141. (3) Pigford, R. L.; Baker, B., III; Blum, D. E. Cycling zone adsorption, a new separation process. Ind. Eng. Chem. Fundam. 1969, 8, 848. (4) Broughton, D.; Gerhold, C. Continuous sorption process employing fixed bed of sorbent and moving inlets and outlets. U.S. Patent No. 2,985,589, 1961. (5) Douglas, J. M. Process Dynamics and Control; Prentice-Hall: Englewood Cliffs, NJ, 1972. (6) Silveston, P. L. Composition Modulation of Catalytic Reactors; Gordon and Breach Science Publishers: Amsterdam, 1998. (7) Silveston, P. L.; Hudgins, R. R. Periodic temperature forcing of catalytic reactions. Chem. Eng. Sci. 2004, 59, 4043-4053. (8) Silveston, P. L.; Hudgins, R. R. Periodic pressure forcing of catalytic reactions. Chem. Eng. Sci. 2004, 59, 4055-4064. (9) Silveston, P. L.; Hanika, J. Challenges for the periodic operation of trickle-bed catalytic reactors. Chem. Eng. Sci. 2002, 57, 3373-3385.

(10) Turco, F.; Hudgins, R. R.; Silveston, P. L. et al. Investigation of periodic operation of a trickle-bed reactor. Chem. Eng. Sci. 2001, 56, 14291434. (11) Tukac, V.; Hanika, J.; Chyba, V. Periodic state of wet oxidation in trickle-bed reactor. Catal. Today 2003, 79-80, 427-431. (12) Kagramanov, Z. G.; Sazonov, A. B.; Magomedbekov, EÄ . P. Simulation of the periodic operating regime of rectifying columns for separating isotopic mixtures of hydrogen. At. Energy (N.Y., NY, U.S.) 2003, 94, 196-199. (13) Bausa, J.; Tsatsaronis, G. Reducing the energy demand of continuous distillation processes by optimal controlled forced periodic operation. Comput. Chem. Eng. 2001, 25, 359-370. (14) Croft, D. T.; LeVan, M. D. Periodic states of adsorption cyclessI. Direct determination and stability. Chem. Eng. Sci. 1994, 11, 1821-1829. (15) Croft, D. T.; LeVan, M. D. Periodic states of adsorption cycless II. Solution spaces and multiplicity. Chem. Eng. Sci. 1994, 11, 1831-1841. (16) Ding, Y.; LeVan, M. D. Periodic states of adsorption cyclessIII. Convergence acceleration for direct determination. Chem. Eng. Sci. 2001, 56, 5217-5230. (17) Ding, Y.; Croft, D. T.; LeVan, M. D. Periodic states of adsorption cyclessIV. Direct optimization. Chem. Eng. Sci. 2002, 57, 4521-4531. (18) Weiner, D. D.; Spina, J. F. Sinusoidal Analysis and Modeling of Weakly Nonlinear Circuits; Van Nostrand Reinhold: New York, 1980. (19) Petkovska, M.; Do, D. D. Nonlinear frequency response of adsorption systems: Isothermal batch and continuous flow adsorber. Chem. Eng. Sci. 1998, 53, 3081-3097. (20) Petkovska, M. Nonlinear Frequency Response of Isothermal Adsorption Controlled by Pore-Surface Diffusion. Bull. Chem. Technol. Macedonia 1999, 18, 149-160. (21) Petkovska, M. Nonlinear frequency response of nonisothermal adsorption controlled by micropore diffusion with variable diffusivity. J. Serb. Chem. Soc. 2000, 65, 939-961. (22) Petkovska, M.; Do, D. D. Use of Higher Order FRFs for Identification of Nonlinear Adsorption Kinetics: Single Mechanisms under Isothermal Conditions. Nonlinear Dyn. 2000, 21, 353-376. (23) Petkovska, M. Nonlinear Frequency Response of Nonisothermal Adsorption Systems. Nonlinear Dyn. 2001, 26, 351-370. (24) Petkovska, M.; Petkovska, L. T. Use of Nonlinear Frequency Response for Discriminating Adsorption Kinetics Mechanisms Resulting with Bimodal Characteristic Functions. Adsorption 2003, 9, 133-142. (25) Guiochon, G.; Shirazi, S. G.; Katti, A. M. Fundamentals of PreparatiVe and Nonlinear Chromatography; Academic Press: London, 1994. (26) Petkovska, M.; Seidel-Morgenstern, A. Nonlinear frequency response of a chromatographic column. Part I: Application to estimation of adsorption isotherms with inflection points. Chem. Eng. Commun. 2005, 192, 1300-1333. (27) Seidel-Morgenstern, A.; Guiochon, G. Thermodynamics of the adsorption of Tro¨ger’s base enantiomers from ethanol on cellulose triacetate. J. Chromatogr. 1993, 631, 37-47. (28) Seidel-Morgenstern, A.; Blumer, C.; Kneip, H. Efficient design of the SMB process based on perturbation method to measure adsorption isotherms and on a rapid solution of the dispersion model. In Fundamentals of Adsorption; Meunier, F., Ed.; Elsevier: Paris, 1998; pp 449-454.

ReceiVed for reView May 19, 2005 ReVised manuscript receiVed October 14, 2005 Accepted October 18, 2005 IE0505965

Fast Estimation of Quasi-Steady States of Cyclic Nonlinear Processes

Recommend Documents