Review Of Binary Distillation - McCabe-Thiele By Computer

Charles R. Koppany, Retired Chemical Engineer

File Name: MTTEXT.DOC
For many years the McCabe-Thiele graphical method (1) has been the accepted procedure for determining
the number of theoretical trays in a binary distillation when constant molal overflow prevails throughout the
column. For close boiling separations such as in an ethane-ethylene or propane-propylene splitter, constant
molal overflow is very closely attained. However, in the rectifying sections of these columns, relative
volatilities are close to one. As a result, the equilibrium curve and upper operating line are very close
together, especially in the region of the overhead product composition xD. This situation makes it very
difficult to graphically step-off theoretical stages, which can be great in number in this region.
The difficulties described above can be circumvented if one is able to solve the McCabe-Thiele method
analytically or by numerical methods on the computer or on a PC. Here a BASIC computer program is
presented for doing just that. This program should be especially useful, for example, to a process engineer,
working on an ethylene plant proposal, who needs to make a quick size (height) estimate for a C2 or C3
splitter column. Use of the program is illustrated for a specific example involving a C 3 splitter.
Finally, a brief discussion is presented for estimating the effect of errors in relative volatility on the
calculation of the number of theoretical stages.
Smoker and Underwood equations. The Smoker (2) and Underwood (3-6) equations are both analytical
techniques for determining the number of theoretical stages in binary distillations operating at conditions of
nearly constant relative volatility and constant molal overflow. However, these procedures do not provide
detailed tray-by-tray liquid and vapor composition profiles that are normally generated by the McCabeThiele (MT).
Figure 1 gives the configuration and nomenclature used for a typical distillation column with a total
condenser and reboiler.

-2Working relationships. First of all, the working relationships used in developing the program will be
briefly described. Attention here will be restricted to single feed, two product binary distillations.
Equilibrium curve. The vapor-liquid equilibrium curve is represented by Eqs. 1 and 2.
= Ax2 + Bx + C

(1)

y = x/{1 + ( - 1)x}

(2)

Here x and y are the mole fractions of the light or lower-boiling component in the liquid and vapor
respectively. is the relative volatility of the light component with respect to the heavier component
(Klt/Khvy). The constants A, B and C are fit to a set of experimental isobaric VLE data, which hopefully
ideally cover the full concentration range for the binary system to be distilled. If VLE data are not
available, then values for the fit will have to be generated by some other means such as from a
generalized correlation which is believed to give reasonably accurate predictions.
Minimum reflux ratio. The minimum reflux ratio (L/D)min is determined by first finding the point were the
q-line intersects the equilibrium curve. This is Point a shown on Fig. 2. The q-line is expressed as:
y = {q/(q - 1)}x + xF/(1 - q)

(3)

q is defined as the moles of liquid introduced to the stripping section resulting from the introduction of each
mole of feed. It basically characterizes the phase condition (mole fraction liquid) of the feed as it enters the
column. Below are summarized the numerical limits for q depending upon various phase conditions of the
entering feed:
Phase Condition

Cold feed
Bubble pt. Fee
Two phase feed

>1
=1
0<q<1

-3Dew pt. Feed

=0

Superheated vapor

<0

Eqs. 1-3 are combined to eliminate y and the resulting expression solved by trial-and-error substitution to
give x/, the liquid mole fraction at the point of intersection of the q-line and the equilibrium curve. y/ is then
calculated from either Eq. 2 or 3. The minimum reflux ratio is then given by
RM = (L/D)min = (xD - y/)/(y/ - x/)

(4)

Upper operating line. The equation of the upper operating line for constant molal overflow is
yn+1 = {R/(R + 1)}xn + xD/(R + 1)

(5)

It is shown as the line segment b-c on Fig. 2. Here R is the actual operating reflux ratio given by some
specified percentage over the minimum, RM. In order to determine the intersection of the upper operating
line and the q-line (Point b), it is necessary to equate Eqs. 3 and 5. The resulting expression is solved
directly for x//, the liquid composition at this point of intersection. y// is then easily determined by either Eq.
3 or 5.
Lower operating line. The lower operating line is now very easily established. It is simply the segment bd with slope given by
Slope = (y// - xB)/(x// - xB)

(6)

Intercept = xB(1 - Slope)

(7)

and intercept by

The program. The computer program, employing all of the above relationships for performing the
McCabe-Thiele graphical construction numerically, is listed in Table 1. The program is written in QBASIC
and has been adapted to run on a PC computer (MS-DOS).
Lines 30-50 request inputs for the coefficients A, B and C of Eq. 1, which along with Eq. 2 establishes the
equilibrium curve. Lines 100-150 require inputs for the product and feed compositions expressed as mole
fraction of the light or lower-boiling component. Lines 160-200 require specification of q, the feed phase
condition, and the ratio of operating to minimum reflux ratio desired. In Line 180 the value of q

-4is set equal to 0.999999 if it was originally inputted as 1.0 (saturated liquid feed). This step avoids the
occurrence of a singularity in Eq. 3.
The necessary parameters generated by Eqs. 1-7 are handled via Lines 210-420. The remainder of the
program is devoted to the stage-by-stage computations based on the previously established equilibrium
curve and the operating lines. Starting with xB, the bottoms liquid composition exiting the reboiler, the
vapor composition yB, exiting the reboiler in equilibrium with xB, is calculated. Obviously, this point (xB,
yB) lies on the equilibrium curve as shown on Fig. 3. From the value of yB the liquid composition from the
bottom theoretical tray is calculated by the equation of the lower operating line. Graphically these
calculated points form the bottom triangle on Fig. 3. This is the first theoretical stage in the calculation
sequence and represents the reboiler. In this way, the program works its way up the column following the
lower operating line and equilibrium curve until the equilibrium vapor composition first exceeds y// (Point
b). At this point the program is ready to switch over to the upper operating line. The number of theoretical
stages (triangles) recorded to this point correspond to the exact number required for the stripping section,
including the feed tray and reboiler. As a result, the feed tray is easily located. It is basically represented
by the triangle with one corner on the lower operating line and the other on the upper operating line.
Finally, the program proceeds to do stage-by-stage computations up the rectifying section by following the
upper operating line. When the liquid composition becomes greater than or equal to x D, the overhead liquid
composition, then the program stops and notes the last stage computed as the top rectifying tray. If the
distillate is a dewpoint vapor, this top tray represents a partial condenser.
The program runs very fast. The speed at which the results are printed out is really the controlling time
constraint.
C3 Splitter example. Let us consider a feed to a C3 splitter containing 0.6954 mole fraction propylene with
the balance being propane. It enters the column as a subcooled liquid with q = 1.2191.

-5It is to be distilled into an overhead liquid product (total condenser) and a bottoms liquid product whose
compositions are 0.9874 and 0.01159 mole fraction propylene respectively. We are asked to determine the
total number of theoretical trays and tray-by-tray composition profiles if the column is to be operated at a
mean pressure of 260 psia and at a reflux ratio (L/D) that is 23.32 % over the minimum.
In 1980 Howat and Swift (7) developed a correlation of propylene-propane VLE data based on the Wilson
liquid phase activity coefficient model for the purpose of determining optimum fractionator operating
pressures in the manufacture of polymer-grade propylene.
The correlation was developed using a rather broad base of all known published VLE data (8-15) for the
propane-propylene system. Employing a database of 216 data points, the authors reported an RMS error of
0.45 % on total pressure. Figure 4 represents a plot of their smoothed values versus pressure from 70 to
130 degrees F over the full range of liquid compositions.
For the C3 splitter application here, values of were read from Fig. 4 along the 260 psia isobar , and then
the coefficients to Eq. 1 were fitted by least squares with a standard deviation, , of 0.01 %.
The tray-by-tray calculations generated from the BASIC program are presented in Table 2. The first
portion of the output summarizes all of the input data. The coefficients to Eq. 1 are printed first followed
then by the minimum and actual reflux ratios and intersection of the operating lines and the q-line. The
remaining output lists stage-by-stage liquid and vapor mole fractions of propylene and the relative
volatilities of propylene with respect to propane. For the sake of brevity, results are printed for only a few
trays around the reboiler, feed tray, and top tray. There are 60 theoretical stripping trays. Tray 1 is the
reboiler and tray 60 is the feed tray. Note here that the vapor from Tray 60, y60=.7205 mole fraction
propylene, just exceeds y//= 0.7177, thus confirming the location of the feed tray. There are 159-61+1 = 99
rectifying trays.
Tray 159 is the top tray because a total condenser is being considered here. The overhead product purity
spec of 0.9874 mole fraction propylene is met because the vapor composition from Tray 159 is 0.9879 mole
fraction propylene.

-6It is interesting to note here that, in the region near the top of the rectifier, the change in composition of
either vapor or liquid is quite gradual- approximately 0.0011 mole fraction/ tray.
This behavior is indicative of how difficult it would be to perform the graphical version of the McCabeThiele method in this section of the column.
Both the Smoker and Underwood correlations were applied to this C3 splitter example. Average relative
volatilities for the rectifying and stripping sections were based on the Howat and Swift VLE data and were
subsequently substituted into the Smoker and Underwood equations to determine the number of theoretical
trays in each of these sections. These values were then compared against the results from the MT program.
The results are quite good and briefly summarized below:

Method

NR

NS

Total

MT Program

99

60

159

Smoker

104

57

161

Underwood

99

56

155

Effect of VLE uncertainties. Following an analysis similar to that conducted by Kister (16) for the effect
of VLE uncertainty on tray efficiency, a brief discussion is presented here for estimating the effect of errors
in relative volatility on the calculation of theoretical trays.
This analysis is generally performed using Fenskes equation for total reflux (minimum theoretical trays).
NM = Ln S/ Ln

(8)

where S, the separation factor, is given by

S = {xD(1 - xB)}/{xB(1 - xD)}

(9)

Attention will be given here to a fixed separation ( S = a constant). As a result, differentiation of Eq. 8
gives
d( NM Ln ) = d Ln S = 0

-7or

(dNM) Ln + NM d Ln = 0

(10)

Upon rearrangement, Eq. 10 becomes

d NM/ NM = - d( Ln )/ Ln

(11)

If the changes or perturbations are small, then to a good approximation

NM/ NM = - ( Ln )/ Ln
or

NM/ NM = -(Ln pred - Ln actual)/ Ln actual

(12)

If it is further assumed that the effect of errors in on NM are the same as on the number of theoretical
trays N at operating reflux, then
N/ N = NM/NM

(13)

Suppose the true or actual value of is 1.1, but that a given VLE correlation predicts high by 5 % such
that pred = (1.05) ( 1.1) = 1.155
Then by Eq. 12 or 13
N/ N = -(Ln 1.155 - Ln 1.1)/ Ln 1.1 = - Ln(1.155/1.1)/Ln (1.1) = -0.512
This means that the number of theoretical trays would be underpredicted by more than 50 %.
Suppose now that the correlation underpredicts by 5 % or that pred = (0.95)(1.1) = 1.045,
then Eq. 12 predicts that Ntheo will be overpredicted by 53.8 %. At lower values of actual, the effect is even
more pronounced.
The brief table below summarizes the effect of a 5 % perturbation in on Ntheo for a range of actual values
typically found in a C3 splitter.
Effect on Ntheo if perturbed by
actual
1.1

+ 5 % error in
-51.2

- 5 % error in
+ 53.8

1.2
1.5

- 26.8
- 12.0

+ 28.1
+ 12.7

-8NOMENCLATURE
A, B and C

Coefficients to Eq. 1

Overhead distillate product rate, moles/hr

Liquid reflux rate from the condenser, moles/hr

Ntheo or N

Total number of theoretical trays at the operating reflux ratio, R

NM

Minimum number of theoretical trays at total reflux

NR

Number of theoretical trays in the rectifying section

NS

Number of theoretical trays in the stripping section

Moles liquid introduced to stripping section resulting from introduction of each mole of
feed

Operating reflux ratio at some percentage over the minimum, L/D

RM

Minimum reflux ratio (infinite stages) or (L/D)min

The separation factor, defined by Eq. 9

Mole fraction of lower boiling or light component in the liquid

x/

Liquid mole fraction of lower boiling or light component at the intersection of

the q-line and equilibrium curve

x//

Liquid mole fraction of lower boiling or light component at the point of intersection of
the operating lines and the q-line

xB

Mole fraction of lower boiling or light component in the bottoms product

xD

Mole fraction of lower boiling or light component in the distillate product

xF

Mole fraction of lower boiling or light component in the feed

Mole fraction of lower boiling or light component in the vapor

yB

Mole fraction of lower boiling or light component in vapor exiting the reboiler,
entering the bottom theoretical tray

-9NOMENCLATURE Continuedy/

Vapor mole fraction of lower boiling or light component at the intersection of

of the q-line and the equilibrium curve

y//

Vapor mole fraction of the lower boiling or light component at the point of intersection
of the operating lines and the q-line

Relative volatility of lower boiling component with respect to the higher boiling
component or equilibrium ratio of lower boiler/equilibrium ratio of higher boiler,
Klt/Khvy

Standard deviation, percent

-10LITERATURE CITED
1.

Mc Cabe, W.L. and E.W. Thiele, Ind. Eng. Chem., 17, No. 6, p. 605, 1925.

2.

Smoker, E.H., Trans. Am. Inst. Chem. Engrs., 34, No. 2, p. 165, 1938.

3.

Underwood, A.J.V., J. Inst. Petrol., 30, p. 225, 1944

4.

IBID., ---------------------------------, 31, p. 111, 1945.

5.

IBID.,----------------------------------, 32, p. 598, 1946.

6.

IBID., Chem. Eng. Prog., 44, p. 603, 1948.

7.

Howat, C.S. and G.W. Swift, Ind. Eng. Chem. Process Des. Dev., 19, p. 318, 1980.

8.

Hakuta, T., Nagahama, K. and M. Hirata, Bull. Jpn. Pet. Inst., 11, p. 10, 1969.

9.

Hanson, G.H., Hogan, R.J., Nelson, W.T., and M.R. Cines., Ind. Eng. Chem., 44, No. 3, p. 604, 1952.

10. Hirata. M., and T. Hakuta., Mem. Tokyo Metropol. Univ., Fac. Technol., 18, p.113, 1968.
11. Hirata, M., Hakuta, T. and T. Onada, Int. Chem. Eng., 8, No. 1, p. 175, 1968.
12. Laurance, D.R. and G.W. Swift, J. Chem. Eng. Data, 17, p. 333, 1972.
13. Manley, D.B. and G.W. Swift, J. Chem. Eng. Data, 16, p. 301, 1971.
14. Mann, A.N., Pardee, W.A., and R.W. Smith, J. Chem. Eng. Data, 8, p. 499, 1963.
15. Reamer, H.H. and B.H. Sage, Ind. Eng. Chem., 43, No. 7, p. 1628, 1951.
16. Kister, H.Z., Distillation Design, Mc Graw-Hill Inc., New York, p. 380, 1992.

Table 1 (Page 11)

Listing of the McCabe-Thiele Tray-By-Tray Program (QBASIC)

10 LPRINT "BINARY DISTILLATION COLUMN DESIGN-SINGLE FEED"
20 LPRINT
30 INPUT "A="; A
40 INPUT "B="; B
50 INPUT "C="; C
60 LPRINT "A", "B", "C"
70 LPRINT A, B, C
80 LPRINT
97 LPRINT
100 INPUT "DISTILLATE COMPOSITION, XD="; XD
110 LPRINT "DISTILLATE COMPOSITION, XD="; XD
120 INPUT "BOTTOMS COMPOSITION, XB="; XB
130 LPRINT "BOTTOMS COMPOSITION, XB="; XB
140 INPUT "FEED COMPOSITION, XF="; XF
150 LPRINT "FEED COMPOSITION, XF="; XF
160 INPUT "FEED THERMAL CONDITION, Q="; Q
170 LPRINT "FEED THERMAL CONDITION, Q="; Q
180 IF Q = 1 THEN Q = .999999
190 INPUT "R/RMIN="; R0
200 LPRINT "R/RMIN="; R0
210 X0 = .5
220 A0 = A * X0 * X0 + B * X0 + C
230 Y0 = (A0 * X0) / (1 + (A0 - 1) * X0)
240 X1 = (Y0 - XF / (1 - Q)) * (Q - 1) / Q
250 E = ABS(X1 - X0)
260 IF E < .0001 THEN 290
270 X0 = X1
280 GOTO 220
290 RM = (XD - Y0) / (Y0 - X0)
300 LPRINT
310 LPRINT "RMIN="; RM
320 R = R0 * RM
330 LPRINT
340 LPRINT "RACTUAL="; R
350 M1 = R / (R + 1)
360 I1 = XD / (R + 1)
370 X2 = (XF / (1 - Q) - I1) / (M1 - Q / (Q - 1))
380 Y2 = M1 * X2 + I1
390 LPRINT
400 LPRINT "INTERSECTION OF O.L.'S AND Q-LINE, X,Y="; X2, Y2
410 M2 = (Y2 - XB) / (X2 - XB)
420 I2 = XB * (1 - M2)
430 N = 1
440 AL = A * XB * XB + B * XB + C
450 YB = (AL * XB) / (1 + (AL - 1) * XB)
460 LPRINT
470 LPRINT "STAGE", "X,LIQUID", "Y,VAPOR", "ALPHA", "TEMP,F"
480 LPRINT N, XB, YB, AL
490 Y = YB
500 X = (Y - I2) / M2

520 AL = A * X * X + B * X + C
530 Y = (AL * X) / (1 + (AL - 1) * X)
540 N = N + 1
550 LPRINT N, X, Y, AL, T
555 IF Y > Y2 THEN 570
560 GOTO 500
570 LPRINT "TOP OF STRIPPING SECTION"
580 LPRINT
590 LPRINT
600 X = (Y - I1) / M1
610 IF X > XD THEN 670
620 AL = A * X * X + B * X + C
630 Y = (AL * X) / (1 + (AL - 1) * X)
640 N = N + 1
650 LPRINT N, X, Y, AL, T
660 GOTO 600
670 LPRINT "TOP OF RECTIFYING SECTION"
999 END
, "TEMP,F"
480 LPRINT N, XB, YB, AL
490 Y = YB
500 X = (Y - I2) /

-12Table 2. Program output for C3 splitter example.

BINARY DISTILLATION COLUMN DESIGN- ANALYTICAL MCCABE-THIELE

A = 9.057862E-04

B = -7.776908E-02

C = 1.167902

DISTILLATE COMPOSITION, XD = 0.9874

BOTTOMS COMPOSITION, XB = 0.01159
FEED COMPOSITION, XF = 0.6954
FEED THERMAL CONDITION, Q = 1.2191
R/RMIN = 1.2332
RMIN = 11.97034
RACTUAL = 14.76182
INTERSECTION OF O. L. S AND Q-LINE, X, Y = 0.6994033

0.7176751

STAGE

X-LIQUID

Y-VAPOR

ALPHA

1
2
3

0.01159
0.01345
0.01555

0.01350
0.01566
0.01810

1.167
1.167
1.167

0.71216
0.71653
0.72054

1.115
1.114
1.114

0.72447
0.72842
0.73239

1.114
1.113
1.113

0.98580
0.98688
0.98794

1.092
1.092
1.092 TOP TRAY

58
0.68939
59
0.69403
60
0.69829
TOP OF STRIPPING SECTION
61
62
63

0.70246
0.70666
0.71088

157
0.98451
158
0.98569
159
0.98685
TOP OF RECTIFYING SECTION

REBOILER

FEED TRAY

-13Future Paper In a future article entitled Review of Binary Distillation McCabe-Thiele by Computer
(Part 2) the analytical procedure (BASIC program) for performing the McCabe-Thiele method will be
extended to handle a two-feed system. In a addition, the capability to calculate theoretical tray
temperatures will also be incorporated.

AUTHORS BACKGROUND

Dr. Charles R. Koppany is a retired chemical engineer formerly employed by C F Braun & Co/ Brown &
Root, Inc. from 1965 to 1994. While at Braun he served in both the Research and Process Engineering
departments. Dr. Koppany has also done part-time teaching in the Chemical Engineering Departments at
Cal Poly University Pomona and the University of Southern California. He holds B.S., M.S. and PhD
degrees in Chemical Engineering from the University of Southern California and is a registered professional
engineer (Chemical) in the state of California.

[email protected]