Computing all the roots of ∑iaiλ−ci=1 equations by a matrix-based approach

Computing all the roots of ∑iaiλ−ci=1 equations by a matrix-based approach

chemical engineering research and design 9 0 ( 2 0 1 2 ) 1697–1700 Contents lists available at SciVerse ScienceDirect Chemical Engineering Research ...

133KB Sizes 5 Downloads 19 Views

chemical engineering research and design 9 0 ( 2 0 1 2 ) 1697–1700

Contents lists available at SciVerse ScienceDirect

Chemical Engineering Research and Design journal homepage: www.elsevier.com/locate/cherd

Short communication

Computing all the roots of matrix-based approach

 a i

i

−ci

= 1 equations by a

Miguel Vacahern a,∗ , Rosendo Monroy-Loperena b a b

División de Ingeniería Química y Bioquímica, Tecnológico de Estudios Superiores de Ecatepec, Av. Tecnológico S/N, 55210 México, México ROMON, Paseo de los Pirules 124, 04250 Distrito Federal, México

a b s t r a c t A useful procedure for finding all the roots of the type of equation indicated in the title, typically used in vapor–liquid separations, is presented. The proposed procedure is the result of a matrix-based analysis to generate a Standard Eigenvalue Problem. The eigenvalues of the resulting matrix system are exactly the roots of the equation indicated in the title. Numerical experiments are presented to show the goodness of the proposed procedure. © 2012 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved. Keywords: Matrix-based Solution; Eigenvalue Problem; Multicomponent distillation; Minimum reflux; Underwood’s equations; Flash vaporization

1.

Introduction

Models for vapor–liquid calculations used extensively in computer-aided process simulation programs employ functions of the type

f () =

N  ai i=1

 − ci

−1=0

(1)

where ai , ci (i = 1, . . ., N) are constants specific for each problem and  is the unknown of the function f(). Examples of this type of calculations are: vapor–liquid flash calculations or short-cut distillation calculations using Underwood’s equations (see for instance Seader and Henley, 2005). Function (1) has a pole at  = ci , that is, a vertical asymptote located at  = ci . Note that N − 1 roots of (1) are located between each pair of ci and one more root can be found in the right side of the pole of greatest value if the ci are arranged so that ci < cj whenever i < j. Since the locations of the poles are known, it is assumed that the seek region is well defined and without problems. However, because the function has discontinuities, coupled with the existence of an inflection point between the poles, the



use of traditional methods such as the Newton’s method frequently diverges or converges to an undesired root or requires an excessive number of iterations. To avoid these convergence problems, approaches such as those of Ripps (1968), that used a combination of half interval and Newton–Raphson; Gritton et al. (2001), that used a homotopic continuation method; Billingsley (2002) and Leibovici and Nichita (2010), that approximated the function, and Koc¸ak (2011), that rearranged the equation and applied a fourth order iterative method, to mention a few, have been proposed. Our main motivation stems to find an efficient, robust, and reliable method to calculate all the roots of Eq. (1). This goal was accomplished by a matrix-based analysis of Eq. (1), as explained in the next section. In the context of vapor–liquid separations in distillation columns, several roots are necessary to calculate by means of Underwood’s equations the minimum reflux and the distribution of components between the column products when the nonkey components are distributed (Seader and Henley, 2005). Also, several roots of Underwood’s equations are needed for direct calculation of pinch points for the rectification or the stripping profiles and to construct distillation spaces (Koehler et al., 1995).

Corresponding author. Tel.: +52 5 5 5000 2324; fax: +52 5 5 5000 2300. E-mail address: [email protected] (M. Vacahern). Received 20 October 2011; Received in revised form 21 March 2012; Accepted 22 March 2012 0263-8762/$ – see front matter © 2012 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cherd.2012.03.013

1698

chemical engineering research and design 9 0 ( 2 0 1 2 ) 1697–1700

The determinant of a rank-one update matrix satisfies (see for instance Meyer, 2000)

Nomenclature a A B c c d D E f I K

constant of Eq. (1) square matrix defined by Eq. (11) square matrix defined by Eq. (12) constant of Eq. (1) column vector column vector square matrix defined by Eq. (6) square matrix defined by Eq. (8) function identity matrix vapor–liquid equilibrium coefficient of a component thermal condition of the feed to the distillation column molar reflux ratio in the distillation column molar reboil ratio in the distillation column eigenvector of Eqs. (14) and (16) mole fraction of a component mole fraction of a component in the feed

q R S x x z

det(D + cdT ) = det(D)(1 + dT D−1 c) Let us define a matrix E as

⎡ ⎢ ⎢ E=⎢ ⎢ ⎣

1 + 1 − ˇ1

1

···

1

1

1 + 2 − ˇ2

···

1

.. .

.. .

.. .

.. .

1

1

(7)

cT = ( 1

1

... 1) (8)

Applying Eq. (5) to Eq. (7) we obtain



N 

det(E) = det(D)(1 + cT D−1 c) =

(j − ˇj ) 1 +

j=1

N  i=1

 1 j − ˇj (9)

then Eq. (3) can be written as

⎡ ⎢ ⎢ det ⎢ ⎢ ⎣

1 + 1 − ˇ1

1

···

1

1

1 + 2 − ˇ2

···

1

.. .

.. .

.. .

.. .

1

1

N 

(j − ˇj ) and after rearranging and

⎤ ⎥ ⎥ ⎥=0 ⎥ ⎦

(10)

· · · 1 + N − ˇN

The main feature of the above representation is that with simple matrix operations we can recast Eq. (2) without troublesome algebraic manipulations. Moreover, if we define the following matrices:



(2)

By multiplying by

· · · 1 + N − ˇN

D = diag(1 − ˇ1 , 2 − ˇ2 , . . . , N − ˇN )

Let us define  i = ci /ai and ˇi = a−1 for i = 1, . . ., N, then Eq. (1) i can be posed as:

i=1

(6)

where

Matrix-based analysis

1 = −1 i − ˇi

⎥ ⎥ ⎥ ⎥ ⎦

E = D + ccT

Subscripts D distillate product of a distillation column i, j integer indices W bottom product of a distillation column

N 



Matrix E can be expressed as a rank-one update matrix

Greek letters ˛ relative volatility of a component ˇ constant of Eq. (2)  constant of Eq. (2) ı constant of Eq. (4)  root of Underwood’s equations  root of Eqs. (1) and (2), eigenvalue vaporization fraction of a vapor–liquid flash

2.

(5)

⎢ ⎢ A=⎢ ⎢ ⎣

1 + 1

1

···

1

1

1 + 2

···

1

.. .

.. .

.. .

.. .

1

1

j=1

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

(11)

· · · 1 + N

equating to zero we obtain N 



(j − ˇj ) 1 +

j=1

N  i=1

and

 1 i − ˇi

=0

⎡ (3)

The left side of Eq. (3) can be written as a determinant using a rank-one update matrix as follows. If D and E are nonsingular matrices of dimension N × N and c and d are N × 1 column vectors, the rank-one update matrix E is E = D + cdT

(4)

ˇ1

⎢0 ⎢ B=⎢ ⎢ .. ⎣ .

0

0

···

0

ˇ2

···

0

.. .

.. .

0



⎥ ⎥ ⎥ .. ⎥ . ⎦

(12)

· · · ˇN

Eq. (10) is represented by det(A − B) = 0

(13)

1699

chemical engineering research and design 9 0 ( 2 0 1 2 ) 1697–1700

Table 1 – Expressions for the solution of Underwood and vaporization flash equations. i

Equation

˛z

ˇi

Interpretation of 

N

i i ˛i −

Feed Underwood’s equation with q = / 1

=1−q

q−1 zi

q−1 ˛i zi

=

=0

− z1

˛i zi

 = − 1

− xR+1

− ˛R+1 x

=

S xi,W

S ˛i xi,W

=

− z1

Ki −1 zi

=

i=1 N

˛z

i i ˛i −

Feed Underwood’s equation with q = 1

i

i=1

N  ˛x

i i,D ˛i −

Rectifying Underwood’ equation

=R+1

i,D

i i,D

i=1 N

˛ x

i i,W ˛i −

Stripping Underwood’s equation

= −S

i=1 N



Flash vaporization equation

zi 1− (1−Ki )

=1

i

i=1

Table 2 – Data and results for the case study 1. i

ai

ci

i

ˇi

1 2 3 4 5 6 7

0.02 0.05 0.10 0.2083. . .3 0.30 0.15 0.10

1.10 1.222. . .222 1.466. . .667 1.833. . .333 2.20 2.75 4.40

55.00 24.44 14.666. . .667 8.80 7.333. . .333 18.333. . .333 44.00

50.00 20.00 10.00 4.80 3.333. . .333 6.666. . .666 10.00

Note that the above system has the form of a General Eigenvalue Problem that can be posed in the following form Ax = Bx

(14)

/ 0 ∈ CN is a where A, B ∈ CN×N ,  ∈ C is an eigenvalue and x = corresponding eigenvector. Equivalently,  is an eigenvalue if det(A − B) = 0. Note that the eigenvalues, , of system (13) are equal to the roots, , of Eq. (2). On the other hand, due to A and B are symmetric, the eigenvalues of the system (9) are real, thus the representation is in agreement with Eq. (2). Then, the determination of the eigenvalues  can be done by applying specific algorithms that exploit the symmetry of matrices A and B (Press et al., 2007). Alternatively, note that B is a diagonal matrix and its inverse is simple. Then system (13) can be recast as det(B−1 A − I) = 0

(15)

having the form of a Standard Eigenvalue Problem, that is B−1 Ax = x

(16)

 1.108326 1.245849 1.515865 1.935004 2.470768 3.073069 4.551674

Table 1 provides appropriate expressions for the solution of Underwood’s equations and one of the vapor–liquid flash equations (Seader and Henley, 2005) by the matrix-based approach of this work.

3.

Example calculations

For illustration purposes, we considered the same problem proposed by Billingsley (2002) and also solved by Leibovici and Nichita (2010) as case study 1, whose parameters ai and ci are shown in Table 2. As case study 2 we considered the solution of Underwood’s feed equation for a mixture of eight components (Seader and Henley, 2005), the relevant input parameters are given in Table 3. The parameters  i and ˇi , are shown in Tables 2 and 3 for both cases. The results obtained from the application of the HQR method (see for instance Press et al., 2007) for both case studies are presented in Tables 2 and 3. The results for case 1 in Table 2 are in perfect agreement with those presented by Billingsley (2002). Table 3 provides the roots of Underwood’s equation for the case study 2. As can be seen, two roots are close to the poles defined by the relative volatility of two components, and then the slope of Eq. (2) is close to

where

⎡ 1 + 1 ⎢ ˇ1 ⎢ 1 ⎢ ⎢ ˇ2 −1 B A=⎢ ⎢ .. ⎢ . ⎣

1 ˇN

1 ˇ1 1 + 2 ˇ2 .. . 1 ˇN

··· ··· .. . ···

1 ˇ1 1 ˇ2 .. . 1 + N ˇN



Table 3 – Data and results for the case study 2.a

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

i

zi

˛i

i

ˇi

1 2 3 4 5 6 7 8

0.0137 0.5113 0.0411 0.0171 0.0262 0.0446 0.3106 0.0354

2.43 1.93 1.00 0.765 0.362 0.164 0.072 0.0362

−9.73723 −0.26090 −3.24574 −7.80117 −5.09160 −2.99103 −0.42949 −3.76836

−4.00709 −0.13518 −3.24574 −10.1976 −14.0652 −18.2380 −5.96516 −104.098

(17)

The solution of Eq. (16) can be accomplished by any conventional solution methods for Standard Eigenvalue Problem systems (see for instance, Saad, 2011).

a

q = 0.8666.

 −6.54635 2.41502 1.04503 0.780140 0.380591 0.184062 0.107159 0.037274

1700

chemical engineering research and design 9 0 ( 2 0 1 2 ) 1697–1700

infinity. Thus, other numerical methods will have difficulties to get convergence.

4.

Conclusions

We propose a practical and reliable method to find all the roots of Eq. (1) or Eq. (2) based on linear algebra approach. In this way the problem of finding all the roots can be seen as a Standard Eigenvalue Problem. This fact enable us to find all the roots of Eq. (1) or Eq. (2) without the need of initial guesses and do not need to make any prevision if some particular root is near a pole of the equations. Moreover, all the roots are found in a single step using Standard Eigenvalue Problem solution techniques.

References Billingsley, D.S., 2002. Iterative solution for

 i

Comput. Chem. Eng. 26 (3), 457–460.

ai −ci

equations.

Gritton, K.S., Seader, J.D., Lin, W.J., 2001. Global homotopy continuation procedures for seeking all roots of a nonlinear equation. Comput. Chem. Eng. 25 (7–8), 1003–1019. Koc¸ak, M.C., 2011. Simple, robust, and fast iterative solution of Underwood’s equation for minimum reflux. Chem. Eng. Res. Des. 89 (2), 197–205. Koehler, J., Poellmann, P., Blass, E., 1995. A review on minimum energy calculations for ideal and nonideal distillations. Ind. Eng. Chem. Res. 34 (4), 1003–1020. Leibovici,  C.F., Nichita, D.V., 2010. Iterative solution for ai −ci

= 1 equations. Chem. Eng. Res. Des. 88 (5–6), 602–605.

i

Meyer, C.D., 2000. Matrix Analysis and Applied Linear Algebra. Society for Industrial and Applied Mathematics. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 2007. Numerical Recipes – The Art of Scientific Computing, third ed. Cambridge University Press. Ripps, D.L., 1968. Minimum reflux figured an easier way. Hydrocarbon Process. (July), 84. Saad, F., 2011. Numerical Methods for Large Eigenvalue Problems, second ed. Society for Industrial and Applied Mathematics. Seader, J.D., Henley, E.J., 2005. Separation Process Principles, second ed. Wiley, New York.