Ind. Eng. Chem. Res. 1993,32, 1767-1771
1767
Comments on the Use of Quadrature in Selecting Pseudocomponents for Multiple-Contact Processes Kraemer D. Luks,' Edward A. Turek, and Tor K. Kragas Amoco Production Company, Box 3385, Tulsa, Oklahoma 74102 Certain complex,multicomponent mixtures can be represented as continuous distributionsof species, typically expressed as a function of one or more parameters (such as carbon number) or physical properties. Rather than use a continuous distribution directly when performing flash calculations, one often uses quadrature to choose a set of discrete components (pseudocomponents) to mimic the distribution. In this paper we critically examine the validity of systematically updating the (quadrature) identification of n pseudocomponents during a multiple-contact process, in order to compensate for the changing compositional nature of the system. We show that updating procedures can produce computationalresults for a multiple-contact process that diverge from those one obtains using n pseudocomponents of fixed identification, where n is very large. As an illustration, the case of a semicontinuous live oil repeatedly contacted with fresh COa is considered.
Introduction It has been shown [e.g., see Cotterman and Prausnitz (1985, 19911, Cotterman et al. (1985a,b, 19861, Behrens and Sandler (1988), and Shibata et al. (198711 that quadrature can be a very efficient way of specifying the pseudocomponents of the C7+portion of a live oil, providing one can accurately represent the distribution of the many C7+ species present as a continuous function of some parameter(s), say, carbon number. For example, given the simple exponential distribution F(I) = ae"'/(e4 - ea) = 0 when I < A, I > B (1) for a fluid portion with ita carbon number I limited to the range A = 6.5 to B = 100.5, one can accurately represent that portion by as few as four pseudocomponenta when performing a single flash calculation. Du and Mansoori (1986,1987) were the first, to our knowledge, to use the simple exponential distribution in phase equilibrium calculations,although Shibata et al. (1987) point out that Vogel et al. (1983)earlier presented experimental evidence supporting the use of the exponential form for C7+ mole fraction distributions. The carbonnumber Ican be related to molecular weight through MW = 14.01+ 2.0 (2) The average carbon number for a given range (A,B) can be expressed as a function of the exponential decay parameter a:
(I)= [ ( a A+ l)e"A - (aB + l)e41/[a(e4
- -
(3) In the limit B a, (I) A + l/a. The higher the molecular weight of the continuous portion, the smaller a becomes, or equivalently, the "flatter" the distribution eq 1becomes. Recently, Matthews et al. (1991) have suggested that, in a multiple-contact process where the continuous portion/ phase composition is changing from stage to stage, one should alter the description of the phase pseudocomponents to reflect that change. In the context of the distribution in eq 1above, alteration would be equivalent to changing the a parameter. It could as well mean a
* Address correspondence to this author at Department of Chemical Engineering, University of Tulsa, 600 S. College Ave., Tulsa, OK 74104-3189.
change in the limit B, but one can argue that a sensible initial choice of B should remain appropriate throughout a multiple-contact process. Matthews et al. carried out their continuous system description using two characterization variables: (1)distillation curve boiling point and (2) specific gravity or Watson characterization factor. They employed a cumulutiue boiling point temperature curve to describe their system, with temperature as the ordinate and cumulative weight fraction as the abscissa. The cumulative weight fraction varies from 0 to 1in every step of the multiplecontact process and therefore permits one to use a Legendre-Gaussquadrature scheme once at the outset of the process to compute the quadrature points and weight factors; the fixed limits (0,l) of the cumulative weight fraction eliminate the need to recompute these for later contactings. The characterization of the pseudocomponenta changeswith successivecontactings,throughchanges in the pseudocomponent properties accompanying the shiftingof the cumulative property curves. In other words, although the locations of the pseudocomponent assignmenta along the cumulatiue property curve and their associated weight fractions conveniently do not change with each contacting, the pseudocomponent properties do. Setting aside any discussion of the relative merits of choosing one set of characterizingparameters over another, we judge that the description employed by Matthews et al. is equivalent in detail to a two-parameterdistribution. In contrast, eq 1 is a one-parameter description. The analysis that follows will focus on certain formal features of Matthews' methodology versus traditional methods using fixed identityfproperty pseudocomponents. The analysiswill be formally,if not quantitatively,independent of description. We will employ the one-parameter exponential function in eq 1in our computations. Ita simplicity, which may properly be considered a quantitative liability, can in fact cause formal differences between Matthews' methodology and more traditional strategies to be numerically magnified. In summary, a recipe analogous to Matthews' can be constructed using the exponential distribution eq 1, even if one must concede certain computational advantages to the use of a cumulative distribution. In the following section, such a computational recipe will be constructed suitable for a distribution function such as presented in eq 1.
0888-588519312632-1767$04.00/0 0 1993 American Chemical Society
1768 Ind. Eng. Chem. Res., Vol. 32,No. 8, 1993
Formalism One can redefine the carbon number variable I in the distribution function as a reduced variable: y=a(I-A)
(4)
Similarly
C = CY(B- A ) (5) represents the upper bound of the reduced variable function. Equation 3 can be reexpressed as (y) = 1 - Ce*/(l- e-‘) = cr((I) -A)
(6) MW is synonymous with (I)at fixed (A$) and directly yields a (and thus C). To select pseudocomponent carbon numbers and mole fractions for a continuous portion of a phase represented by this distribution, one can apply finite Laguerre-Gauss quadrature, which is applicable to integrals of the type S,‘f(z)e+ d~
(7)
The quadrature weighting factors {wi)and the roots {xi) of the pseudocomponents are functions of the reduced parameter C. {Ii} depend directly on {xi]through eq 4. If the continuous portion of a fluid phase has a MW that varies with the number of process steps performed, then {wi]and (Ii)can in turn be varied with the steps through C. Since it is understood that the properties of a pseudocomponentdepend on its carbon number,then they will also vary with the steps. In the limit C m (Le., B Q),given a fixed value of A), the quadrature results {wj]and {xi)will be fixed at all process steps for the continuous portion of a phase. The MW (equivalently the value of (I))of the continuous portion at each process step will determine a value of a at that step, and consequently the {Ij]will vary with each step along with MW. This limiting scenario is mathematically simple, like Matthews’ which employs a cumulative distribution function accompanied by LegendreGauss quadrature. But it is not formally equivalent to the Matthews methodology. Matthews et al. were in essence able to achieve this degree of mathematical simplicity for arbitrary values of C; that is, they were able to impose a finite upper bound B on the carbon numbers of the continuous portion at the outset and calculate their Legendre-Gauss quadrature results once (initially), suitable for all later process steps. For finite B in the context of the exponential distribution in eq 1,one is required to recompute {wi)and (xi) at each step of the process. This is a computational inconvenience, of course. Of importance here, however, is that this computation for a value of C that varies with each step, inconvenience aside, is formally equivalent to that of Matthews et al., as B is maintained aa a finite, fixed parameter for the entire multiple-contact process. We will now compare the computational results of this last recipe,cumbersomethough it may be, to those of more traditional strategies whereby the identity and properties of the pseudocomponentsare initiallyestablishedand then held fixed throughout the multiple-contact process. By this comparison, we will be able to evaluate the ramifications of Matthews’ methodology relative to traditional methodologies within the framework of a common distribution, that being the one in eq 1.
-
-
Example Calculation A live oil, ‘oil A” (described in Table I), with 60.7 mol 5% C7+was chosen for the initial liquid phase of a multiple-
Table I. Molar Composition of Oil A Used in the Example Calculations component mole fraction component mole fraction N2 0.0216 iC4 0.0188 c1 0.1080 nC4 0.0338 COa 0.0187 iC6 0.0215 Ca 0.0520 nC5 0.0244 H2S 0.0024 CS 0.0457 cs 0.0461 C7+ 0.6070
contacting proms. Two situations were examined, wherein the molecular weight of the C7+ portion of oil A was arbitrarily assigned to be either 149or 249. The first value is close to that used by Matthews et al. (1991)for his CS+ mixture. The following multiple-contact processes were computed at 160 OF using fresh C02 as the contacting gas, with the C7+ portion represented by eight (quadrature) pseudocomponents: 1. Fifteen contactings were made of the liquid phase, initially oil A having a C7+ MW of 149,with equimolar amounts of C02 at 1500 psia, after which the equilibrium vapor phase was removed; four different values of the carbon number upper bound B of the C7+ portion were examined, 20.5, 40.5,60.5,and m, respectively. Carbon number lower bound A is fixed at 6.5in all computations. 2. Fifteen contactings were made of the liquid phase, initially oil A having a C7+ MW of 249,with twice the molar amount of C02 at 2500 psia, after which the equilibrium vapor phase was removed; three different values of the carbon number upper bound B of the C7+ portion were examined, 60.5,80.5, and m, respectively. Carbon number lower bound A is fixed at 6.5 in all computations. 3. Fifteen contactings analogous to 1 were performed but with the eight C7+pseudocomponent carbon numbers (Le., properties) determined once, prior to the first contacting, and held constant. 4. Fifteen contactings analogous to 2 were performed but with the eight C7+pseudocomponent carbon numbers (i.e., properties) determined once, prior to the first contacting, and held constant. As the MW of the liquid phase increased with each contacting in 1 and 2,the carbon numbers of the eight C7+ pseudocomponents were readjusted prior to the next contacting by finite Laguerre-Gauss quadrature to be consistent with that MW and the (fixed) carbon number upper bound B. In performing our computations, we empiricallyfound that our multiple-contactprocess results did not vary significantly when we used six or more quadrature components. Therefore, we arbitrarily used eight components for our computations. When applying Laguerre-Gauss quadrature to cases of more than six pseudocomponents,we invoked quadruple precisionwhen generating the roots and weighting factors. We deemed this necessary from comparisons of calculated roots and weighting factors with literature tabulations for the case C m (Abramowitzand Stegun, 1972). All computations were performed using the same Redlich-Kwong equation of statedescription as employed in our earlier paper (Luks et al., 1990). It should be emphasized that the choice of equation of state is not crucial here, only that the same equation of state description be used in comparative computations. Larger values of B as upper carbon number bounds were employed for the computations for MW = 249 than for those for MW = 149. MW = 249 corresponds to a carbon number of about 17.6,versus 10.5 for MW = 149. When C ( B a),the normalized C7+mole fractions are the same for all contactings, and thus quadrature only has to be performed once; however, the properties of the
-
- -
Ind. Eng. Chem. Flea., Vol. 32, No.8,1993 1769
e
460
I
I
L
.-c1)
Sliding (24W60.6)
Discrete (240/60 6)
al
Sliding (14W20.5) Sliding (140/B -. -) Discrete (149/20.6)
Discrete (140/B + m) I
60
0
3
6
9
3
12
3
0
16
pseudocomponents are changing with each contacting as the carbon numbers change. In their work, Matthews et al. performed a differential pressure depletion calculation consisting of four steps at 533.2 K for a mixture of 30 wt 7% CO2 + 70 wt 7% c6+ mixture. Their c6+ mixture had a molecular weight of about 145. The first depletion removed a substantial portion of the pressuring gas, similar to what we experienced with the light ends of the live oil in our computations. It is not straightforward to assess how similar the multiple contact process calculations 1 and 2 above are to Matthew’ oil depletion calculation. Such an assessment, however, is not needed. The computationsperformed herein focus on their methodology of adjusting pseudocomponent properties at each step of a process vis-a-vis more traditional schemes using pseudocomponents with fued properties. The comparison is made in a common framework, the above-listed multiple-contact processes.
I
I
6
9
12
16
Number of Contactings
Number of Contactings Figure 1. Comparieon of the liquid-phase molecular weight a~ a function of the number of contactings for the sliding and discrete methods for two seta of initial CI+ molecular weight and parameter B.
I
BB a function of the numbr of contacting8 for the sliding and discrete meth~fortwo~teofinitialCI+molecular~ightandtheparame~ B 0,
Figure 2. Comparison of the liquid-phase molecular weight
-
2oo
3
50
-
Sliding (149/20.5)
3 I
0
I
I
1
Discussion For convenience, we will hereafter refer to the characterization updating scheme used in calculations 1 and 2 as the sliding scheme, since the carbon numbers of the remaining liquid-phase pseudocomponents in these cases are adjusted (“slide”)to higher carbon numbers with each contacting. This sliding will be most pronounced for the case B m, as the freedom for the carbon numbers to adjust to the changing value of MWl is not restricted by a finite carbon number upper bound B. The characterization scheme where the carbonnumbers (and properties) are fixed initially and held constant throughout the multiple contactings will be referred to as the discrete scheme. Figures1 and2compareMWlfortheslidinganddiscrete computations for the cases of CI+initial MW = 149 for B = 20.5 and m, and C7+ initial MW = 249 for B = 60.5 and m. The computations for MW = 249 with the different values of E diverge with the number of contactings more dramatically than those for MW = 149.(There is a “kink” in all the curves at contacting 1,more noticeable for C7+ initial MW = 249,which is caused by the first contacting preferentially removing nearly a l l the volatile8 of the original live oil.) The same type of property divergence (but not as pronounced) occurs when one computes the CI+mass left in the liquid phase (Figures 3 and 4). The divergences between sliding and discrete results for MW = 149 are very small in all the Figures 1-4 when compared
160
-
Discrete (140; B -. 4 ._
Sliding (140; B -. m)
60-
0
3
6
Q
12
15
Number of Contactings Figurn 4. Comparison of the liquid-pham mam aa a function of the number of contactingn for the sliding and discrete methodn for two seta of initial Cy+ molecular weight and the parameter B =,
-
with MW = 249,especially over the fmt five contactings. The divergences increase not only with the number of contactings but ale0 with the value of the upper carbon number bound B. For an initial C7+ MW = 149, for a particular computational strategy (sliding or discrete), the resulta for B 1 40.5 are virtually identical to those for
1770 Ind. Eng. Chem. Res., Vol. 32,No. 8,1993 0.5
1
Legend
8
0.3
c
0.5 I
Legend
C
Initial (148/20.6)
.-50
0.4
Sliding (16) Discrete (161
8
0.3
-
c
Initial (249/60.6)
/--\
Sliding (16) Discrete (16)
\
-
2
x
s
0.2
0.1
3
0.0 1
2
3
4
5
6
7
8
C7+ Pseudocomponent c7+mole fractions (normalized)for ininitialc,+molecular
pigun, S, weight of 149 and B = 20.5. Shown are the initial values, the values after 16 contactingafor the sliding and discrete calculations,and the values for B a.
-
I
0.5 I
Legend Initial (14W606)
slifing(16) - Discrete (16) Sliding (B -1 +
1
2
3
4
5
6
7
8
C7+ Pseudocomponent Fieurn 6. C7+ mole fractions (normalized)for an initialC7+molecular weight of 149 and B = 60.5. Shown are the initial values, the values after 15 contactinge for the sliding and discrete calculations,and the values for B a.
-
-.
B a at five contactings, and likewise for B 1 60.5 a t 15 contactings. For an initial C7+MW = 249,it appears that one has to impose a bound B > 80.5 to accomplish this. Given the above results, we wish to address two questions: 1. Why do the sliding and discrete results differ? 2. Which result is "correct"? (It is recognized that whereas neither result may be "exact", one may be formally more correct than the other.) Some preliminary insights to the first question can be obtained by examining how the distribution of the normalized C7+ mole fractions shifts as the coatactings proceed. Figures 5 and 6 illustrate these distributions for the initial C7+ MW = 149, for B = 20.5 and 60.5, respectively. Smoothed curves have been conatructed through the eight normalized mole fractions in each case. Interestingly, the initial mole fractions of the pseudocomponents that represent an exponentialdistribution exhibit a maximum, which is at the second or third pseudocomponent. For f i i t e valuesof B, both the discrete and sliding mole fraction distributions shift in a logical manner, with the presence of the heavier pseudocomponents increasing at the expense of the lighter ones. The "Sliding ( B 05)" mole fraction curve is invariant with the number of contactings for the sliding computation, as pointed out earlier.
-
1
2
3
4
5
6
7
8
C7+ Pseudocomponent Figun, 7. C7+mole fractions(normalizsd)for an initial C7+molecular weight of 249 and B = 60.6. Shown are the initial values, the values after 16 contactingafor the sliding and discrete calculations,and the values for B m.
-
The distribution for the sliding case changes form with the number of contactings differently from that for the discrete case. Bothstart with the same initial distribution (shown), with all the pseudocomponents of the sliding case experiencing increasing carbon numbers with the number of contactings. One effect of these carbon number changes is a reduction in the volatility of the lighter C7+ pseudocomponents, causing their normalized mole fractions not to diminish as rapidly as in the discrete case. The sliding pseudocomponent mole fraction distribution thus becomes flatter with the number of contactingswhen compared to the corresponding discrete distribution. Since the pseudocomponents are changing identity (carbonnumber) with each contacting for the sliding case, one cannot strictly compare pseudocomponent masses between contactings with the objective of quantitatively examining a possible conservation-of-mass violation. However, it was noted that the amount of mass in the heaviest pseudocomponents can and will increase for the sliding case when the limit B is finite. For example, this occurred with the sixth through eighth C7+ pseudocomponents for initial C7+MW = 149 when the upper carbon number bound B = 20.5and 40.5,and with the f i i h through eighth C7+ pseudocomponents when B = 60.5. This anomaly, which cannot occur with the discrete case, disappears for the sliding case as B a. In comparing Figures 5 and 6,there is not as much a shift in the normalized mole fraction curve peak with pseudocomponent number for upper carbon number bound B = 60.5 as for B = 20.5. The sliding case distributions for finite B's will be closer to that of the case "Sliding (B m)" as B gets larger, which is logically to be expected. The sliding distribution is flatter after 15 contactings than the discrete distribution for B = 60.5, as before with B = 20.5. Figures 7 and 8,with the C7+ initial MW now equal to 249,are analogous to Figures 5 and 6. The differences in the evolution of the pseudocomponent normalized mole fraction distributions with contactings are more pronounced, which is consistent with the larger property differences seen in Figures 1-4. Again, the distributions for the sliding cases are flatter than those for the discrete cases. The dramatic divergences of results in Figures 1-4 are therefore not surprising. We can now answer the fiist question as to why the discrete and sliding results differ. The mechanics of the discrete computation are straightforward. Given an assumed initial exponentialdistribution of the C7+portion
-
-
Ind. Eng. Chem. Res., Vol. 32, No. 8, 1993 1771 0.5 I
1
are already becoming prominent at five contactings and are more sensitive to the C7+upper carbon number bound
B. Nomenclature
3 U
0.1
t' / I
\
/'
2
1
3
4
\
5
--_ 0
I
7
8
C7+ Pseudocomponent Figure8. C7+mole fractions (normalized)for an initial C7+molecular weight of 249 and E = 80.5. Shown are the initial values, the values after 15 contactinga for the sliding and discrete calculations, and the values for E m. +
of oil A, one gets a distribution of normalized C7+ mole fractions which shifts with the number of contactings. The distribution of these normalized mole fractions is not consistent any longer with an exponential distribution after the first contacting. The MWI information is derived information and not used to adjust later C7+distributions. In contrast, in the sliding property approach, the C7+ distribution is always designated to be exponential for the liquid phase (tobe contacted) and enters each equilibrium contacting coupled to the MWIof the previous contacting. This is an additional constraint on the computation that creates a bias in the computation, one that in the context of multiple contactings causes divergences ( as illustrated in Figures 1-4) between the computed results of the two methodologies. We now address the second question concerning which result is "correct". It is not just the contacting-tocontacting constraint of a specified distribution (herein, exponential) being applied that leads us to believe that the sliding approach is flawed. Intuitively, it is an appealing approach. However, as pointed out earlier, when one uses six C7+pseudocomponents derived by quadrature techniques in the discrete approach for a problem consistingof 15contactings,one gets results that are essentially identical with those that would be obtained using more than six pseudocomponents, for example, using an indefinitely large number of pseudocomponents that could be envisioned to approach infinity in number. Given the premise that the initial exponential distribution correctly represents the extremely large number of C7+ species in a "laboratory" system, then the discrete answer for n L 6 (--) is the correct answer to the multiple-contacting "laboratory" problem. The sliding approach produces an answer that does not asymptotically approachthis discrete answer for n 1 6 (+oJ). The calculation of Matthews et al. (1991)is not a severe test of the methodology that they propose, and therefore any difficulty that may arise from it is not apparent. Their calculation is similar in severity to the computations for the C7+ upper carbon number bound B = 20.5 in Figures 1 and 3 for MW = 149. Very little change in the pseudocomponent carbon numbers takes place with each contacting, and the sliding and discrete computationsare essentially indistinguishable at five contactings. This is not the case for a greater number of contactings and/or for higher C7+ MW's where the deviations in Figures 1-4
A = lower carbon number bound B = upper carbon number bound C = upper bound of the reduced carbon number variable F = distribution function I = carbon number MW = molecular weight w = quadrature weighting factor x = quadrature root y = reduced carbon number variable z = dummy variable cy = exponential distribution decay parameter Subscripts i = species 1 = liquid phase Literature Cited Abramowitz, M., and Stegun, 1. A., Eds. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables; Dept. of Commerce: Washington, DC, 1972; loth printing with corrections, p 923. Behrens, R. A.; Sandler, S. I. The Use of SemicontinuousDescriptions to Model the C7+ Fraction in Equation of State Calculations. SPE Reservoir Eng. 1988,3,1041-1047. Cotterman, R. L.;Prausnitz, J. M. Flash calculations for Continuous or Semicontinuous Mixtures Using an Equation of State. Znd. Eng. Chem. Process Des. Dev. 1986,24,434-443. Cotterman, R. L.;Prausnitz, J. M. Continuous Thermodynamics for Phase-Equilibrium Calculations in Chemical Procese Design. In Kinetic and Thermodynamic Lumping of Multicomponent Mixtures; Astarita, G.,Sandler, S., E&.; Elsevier: Amsterdam, 1991; pp 229-275. Cotterman, R. L.;Bender, R.; Prausnitz, J. M. Phase Equilibria for Mixtures Containing Very Many Components. Development and Application of Continuous Thermodynamics for ChemicalProcess Design. Znd. Eng. Chem. Process Des. Dev. 198Sa,24,194-203. Cotterman, R. L.; Dimetrelis, D.; Prausnitz, J. M. Design of Supercritical-Fluid-ExtractionProcesses Using Continuous Thermodynamics. In Supercritical Fluid Technology; Penninger, J. M. L.,Radosz, M., McHugh, M. A,, Krukonis, V. J., Eds.; Elsevier: Amsterdam, 1985b;pp 107-120. Cotterman, R. L.;Chou, G. F.; Prausnitz, J. M. Comments on "Flash Calculations for Continuous or Semicontinuous Mixtures Using an Equation of State". Znd. Eng. Chem. Process Des. Dev. 1986, 25,840-841. Du, P. C.; Mansoori, G. A. Phase Equilibrium Computational Algorithms of Continuous Mixtures. Fluid Phase Equilib. 1986, 30, 57-64. Du, P. C.; Mansoori, G. A. Phase Equilibrium of Multicomponent Mixtures: Continuous Mixture Gibbs Free Energy Minimization and Phase Rule. Chem. Eng. Commun. 1987,54,139-148. Luke,K. D.; Turek, E. A.; Kragas, T. K. Asymptotic Effecta Using Semicontinuous vis-a-vis Discrete Descriptions in Phase Equilibrium Computations. Znd. Eng. Chem. Res. 1990,29,2101-2106. Matthews, M. A.; Mani, K. C.; Haynes, H. W. Continuous Phase Equilibrium Thermodynamics for Sequential Operations. In Kinetic and Thermodynamic Lumping of Multicomponent Mixtures; Astarita, G., Sandler, S., Eds.; Elsevier: Amsterdam, 1991; pp 307-324. Shibata, S. K.; Sandler, S. I.; Behrens, R. A. Phase Equilibrium Calculationsfor Continuous and SemicontinuousMixturea. Chem. Eng. Sci. 1987,42,1977-1988. Vogel, J.L.;Turek,E. A.; Metcalfe, R. S.; Bergman, D. F. Applications of Equations of State to Calculate Reservoir Fluid Properties. Fluid Phase Equilib. 1983,14,103-116. Received for review November 16, 1992 Accepted April 27,1993
