# Equationtearing Procedures Using The Tridiagonalmatrix Algorithm

As seen earlier, the manual Thiele-Geddes method involves solving the equilibrium-stage equations one at a time. More powerful, flexible, and reliable computer programs are based on the application of sparse matrix methods for solving simultaneously all or at least some of the equations. For cases in which combined column feeds represent mixtures that boil within either a narrow range (typical of many distillation operations) or a wide range (typical of absorbers and strippers) and in which great flexibility of problem specifications is not required, equation-tearing procedures that involve solving simultaneously certain subsets of the equations can be applied. Two such equation-tearing procedures are the bubble-point (BP) method for narrow-boiling mixtures suggested by Friday and Smith (op. cit.) and developed in detail by Wang and Henke [Hydrocarbon Process., 45 (8), 155 (1966)], and the sum-rates (SR) method for wide-boiling mixtures proposed by Sujuta [Hydrocarbon Process., 40(12), 137 (1961)] and further developed by Burningham and Otto [Hydrocarbon Process., 46(10), 163 (1967)]. Both methods are described here. However, the BP method has been largely superseded by the more reliable and efficient simultaneous-correction and inside-out methods, which do, however, incorporate certain features of the BP method. Both the BP and SR methods start with the same primitive equations for the theoretical model of an equilibrium stage as presented next.

Consider a general, continuous-flow, steady-state, multicomponent, multistage separation operation. Assume that phase equilibrium between an exiting vapor phase and a single exiting liquid phase is achieved at each stage, that no chemical reactions occur, and that neither of the exiting phases entrains the other phase. A general schematic representation of such a stagej is shown in Fig. 13-48. Entering stage j is a single- or two-phase feed at molal flow rate Fj, temperature TFj, and pressure PFj and with overall composition in mole fractions zij. Also entering stage j is interstage liquid from adjacent stage j _ 1 above at molal flow rate Lj _ ^ temperature Tj _1, pressure Pj _ ^ and mole fractions %ij _!. Similarly, interstage vapor from adjacent stage j + 1 below enters at molal flow rate Vj+!, Tj+!, Pj+1 and mole fractions yij+^ Heat is transferred from (+) or to (_) stage j at rate Qj to simulate a condenser, reboiler, intercooler, interheater, etc. Equilibrium vapor and liquid phases leave stage j at Tj and Pj and with mole fractions y,j and xij respectively. The vapor may be partially withdrawn from the column as a sidestream at a molal flow rate Wj, with the remainder Vj sent to adjacent stage j _ 1 above. Similarly, exiting liquid may be split into a sidestream at a molal flow rate of Uj, with the remainder Lj sent to adjacent stage j below.

For each stage j, the following 2C + 3 component material-balance (M), phase-equilibrium (E), mole-fraction-summation (S), and energy-balance (H) equations apply, where C is the number of chemical species:

_ 1xi,j _ 1 + Vj+1yi,j+1 + Fjzij _ (Lj + Uj)Xi,j _ (Vj + Wj)yij = 0 (13-68) y,j _ KtjXtj = 0 (13-69)

In general, K values and molal enthalpies in these MESH equations are complex implicit functions of stage temperature, stage pressure, and equilibrium mole fractions:

where vectors Xj and yj refer to all i values of x,j and yi,j for the particular stagej. As shown in Fig. 13-49, a general countercurrent-flow col FIG. 13-48 General equilibrium stage.

umn of N stages can be formed from a collection of equilibrium stages of the type in Fig. 13-48. Note that streams L0, VN+r, W1 and UN are zero and do not appear in Fig. 13-49. Such a column is represented by N(2C + 3) MESH equations in [N(3C + 10) + 1] variables, and the difference or [N(C + 7) + 1] variables must be specified. If the specified variables are the value of N and all values of z,j, Fj, TFj, PFj, Pj, Uj, Wj, and Qj, then the remaining N(2C + 3) unknowns are all values of y,j, %j, Lj, Vj, and Tj. In this case, Eqs. (13-73), (13-74), and (13-77) are nonlinear in the unknowns and the MESH equations can not be solved directly. Even if a different set of variable specifications is made, the MESH equations still remain predominantly nonlinear in the unknowns. For the BP method as applied to distillation, specified variables are those listed except that bottoms rate LN is specified rather than partial reboiler duty QN. This is equivalent by overall material balance to specifying vapor-distillate rate V1 in the case of a partial condenser or liquid-distillate rate U1 in the case of a total condenser. Also, reflux rate L1 is specified rather than condenser duty Q1. For the SR method as applied to absorption and stripping, the specified variables are those listed without exception.

Tridiagonal-Matrix Algorithm Both the BP and the SR equation-tearing methods compute liquid-phase mole fractions in the same way by first developing linear matrix equations in a manner shown by Amundson and Pontinen [Ind. Eng. Chem., 50, 730 (1958)]. Equations (13-69) and (13-68) are combined to eliminate y,,j and y,,j+1 (however, the vector yj still remains implicitly in K,,j):

Next, Eq. (13-68) is summed over the C components and over stages 1 through j and combined with Eqs. (13-70), (13-71), and "Z, z,j -

1.0 = 0 to give a total material balance over stages 1 through j:

By combining Eq. (13-73) with Eq. (13-74), Lj is eliminated to give the following working equations for component material balances: FIG. 13-49 General countercurrent cascade of N stages.

where Aj — Vj + I (Fm - Wm - Um) - V, 2 < j < N (13-76) Bi,j — - Vj + 1 + I (Fm - Wm - Um) - V + Uj + (Vj + Wj)K,j

1 << j << N (13-77) 1 << j << N - 1 (13-78)

The NC equations (13-75) are linearized in terms of the NC unknowns xij by selecting unknowns Vj and Tj as tear variables and using values of vectors Xj and yj from the previous iteration to compute values of K,j for the current iteration. In this manner all values of Aj, B,j, and Ctj can be estimated. Values of Dt j are fixed by feed specifications. Furthermore, the NC equations (13-75) can be partitioned into C sets, one for each component, and solved separately for values of xi, which pertains to all j values of xi,j for the particular species i. Each set of N equations is a special type of sparse matrix equation called a tridiagonal-matrix equation, which has the form shown in Fig. 13-50a for a five-stage example in which, for convenience, the subscript i has been dropped from the coefficients B, C, and D. For this type of sparse matrix equation, we can apply a highly efficient version of the gaussian elimination procedure called the Thomas algorithm, which avoids matrix inversion, eliminates the need to store the zero coefficients in the matrix, almost always avoids buildup of truncation errors, and rarely produces negative values of x,j.

The Thomas algorithm begins by a forward elimination, row by row starting down from the top row (j = 1, the condenser stage), to give the following replacements shown in Fig. 13-50b. For row 1:

where ^ means "is replaced by."

For all subsequent rows:

Pj = Cj'(Bj - AjPj-1), qj = (Dj - Aq- 1)'(Bj - AjPj_ 1),

At the bottom row for component i, xN = qN. The remaining values of Xj for species i are computed recursively by backward substitution:

BP Method for Distillation The bubble-point method for distillation, particularly when the components involved cover a relatively narrow range of volatility, proceeds iteratively by the following steps, where k is the iteration index for the entire distillation column.

1. Specify N and all values of z,j, Fj, TFj, PFj, Pj, Uj, Wj, and Qj, except Q1 and QN.

2. Specify type of condenser. If total (U1 ^ 0), compute LN from overall material balance; if partial (U1 = 0), specify V1 and compute LN from overall material balance.

3. Specify reflux rate L1, assuming no subcooling.

5. Provide initial guesses (k = 0) or values of all tear variables Tj and Vj (j > 2). Temperature guesses are readily obtained by linear interpolation between estimates of top- and bottom-stage temperatures. The bottom-stage temperature is estimated by making a bubble-point-temperature calculation by using an estimate of bottoms composition at the specified bottom-stage pressure. A similar calculation is made at the top stage by using an estimate of distillate composition; otherwise, for a partial condenser, a dew-point temperature calculation is made. An estimate of the vapor-rate profile is readily obtained by assuming constant molal overflow down the column.

6. Set index k = 1 to initiate the first column iteration.

7. Using specified stage pressures, current estimates of stage temperatures, and current estimates of stage vapor- and liquid-phase

 B, c,