ARTICLE pubs.acs.org/IECR
An Object-Oriented Approach for Structure Design of Property Calculation Programs Using Equations of State Lijuan Shi*,† and Xiaodong Liang‡ † ‡
School of Biological and Chemical Engineering, Jiaxing University, Jiaxing, Zhejiang 314001, People’s Republic of China Honeywell Technology Solution (China), 430 Li Bing Road, Pudong New Area, Shanghai 201203, People’s Republic of China ABSTRACT: A new object-oriented approach for designing general property calculation programs using different types of equations of state (EoS) is presented in this work. It helps to ease the programs’ extensibility, reusability, and maintainability and, at the same time, minimize the code redundancy. The approach consists of two parts: (a) the total Helmholtz free energy function is separated into specific terms, based on physical meanings or mathematical expressions, independent of the properties of interest, and the same or similar specific terms used in different equations of state are abstracted to general algorithms; (b) the general relationships between these algorithms and the properties of interest are implemented at the highest level of the program. According to this strategy, a general cubic algorithm can be abstracted for all of the cubic EoS. The general hard-sphere, hard-chain, and association algorithms can be aggregated into different Statistical Associating Fluid Theory (SAFT) family EoS with different dispersion algorithms. Similarly, the general cubic and association algorithms can be aggregated into Cubic Plus Association (CPA) EoS. Taking Perturbed-Chain SAFT (PC-SAFT) EoS as an example, this manuscript demonstrates a clear object-oriented structure for developers, with regard to how to design a complete thermodynamic model package using this approach.
1. INTRODUCTION Equations of state (EoS) play an important role in numerous industries, and they have assumed an expanding role in the study of phase equilibria, volumetric behavior, and thermophysical property calculations for fluids, fluid mixtures, and solids, etc. They are widely used in presenting vaporliquid, liquidliquid, supercritical fluidliquid/solid and other types of phase equilibrium behaviors over wide ranges of temperature and pressure, and mixtures of diverse compounds. Over the past five decades or so, many EoS have been developed to describe the thermodynamic behaviors of fluids. Soave RedlichKwong (SRK)1 and PengRobinson (PR)2 EoS are two of the most popular methods, especially for hydrocarbon systems. They have been also improved and extended to wider applications. Among them, the PengRobinson StryjekVera (PRSV)3 with empirical asymmetric mixing kij formulation and TwuSim Tassone (TST)4,5 with excess Helmholtz free-energy mixing rule are two successful variants. Recently, a trend has emerged where new EoS are developed, based more on theoretical considerations. Generally, the dependent quantity—usually pressure or Helmholtz energy—is assumed to be comprised of a series of terms, each accounting for a particular thermodynamic phenomenon or behavior. Often, the equations include a repulsive reference term that is based on the results of molecular simulations, using some type of simplified potential function (hard sphere, hard chains, or LennardJones, for example). Other terms then account for phenomena such as association and polarity. Once the terms for all desired phenomena have been included, the EoS is established. Perturbed HardSphereChain (PHSC)6,7 and EoS based on Statistical Associating Fluid Theory (SAFT), such as CK-SAFT,8,9 SAFT-VR,10 SAFTLJ,11,12 PC-SAFT,1315 GC-PC-SAFT,16 etc., are all in this category. On the other hand, introductions of theoretical terms into cubic equations of state (CEOS) could be tempting, because CEOS are r 2011 American Chemical Society
easy to solve and have accumulated vast industry applications throughout the years. The Cubic Plus Association1719 (CPA) EoS is one of the most successful models, which has been successfully applied to systems containing alcohols, glycols, water, alkanes, organic acids, and amines. Although many EoS have been developed for different types of systems; many of them share the same general structure or have the same terms. For example, different variants of CEOS share the same general structure with mere differences in parameters and/or mixing rules; CK-SAFT8,9 and PC-SAFT13,14 are different from the formulas merely in the dispersion term, and CPA17 is a CEOS plus an association term, which can be shared with the SAFT family EoS. If each model has its own complete code implementations in the program, it will make the bloated program difficult to extend, reuse, and maintain. Thus, in an ideal situation, one general algorithm would have only one code copy in the entire program, and different EoS would be given the flexibility of having their own alternative implementations for special treatments. In this work, a general object-oriented (OO) structure is proposed to achieve this target, which could be used in large commercial process simulators, as well as standalone thermodynamic calculation tools. Some attempts to design an OO simulator and/or standalone thermodynamic tools have been made,20,21 which focused more on the design of interfaces for different components, while this work totally focuses on the internal thermodynamic package structure. The feasibility of this approach is analyzed in the General Structure section. A detailed property package structure design is described in the Specific Received: July 20, 2010 Accepted: March 30, 2011 Revised: March 16, 2011 Published: April 12, 2011 6404
dx.doi.org/10.1021/ie101547t | Ind. Eng. Chem. Res. 2011, 50, 6404–6412
Industrial & Engineering Chemistry Research
ARTICLE
Package Structure Design section, using PC-SAFT as an example. The schematic structure of the entire program, and the benefits of this approach, are highlighted in the Results and Conclusion sections. Useful expressions of the PC-SAFT EoS are given in the Appendix.
2. GENERAL STRUCTURE From the general thermodynamic point of view, it is easy to develop consistent models and calculations for properties by combining the residual Helmholtz free energy (Ar) and its partial derivatives. This is because the surface, which is composed of the Helmholtz free energy (A), temperature (T), total volume (V), and number of moles (n), can completely describe the thermodynamic systems. The properties of engineering interest that can be obtained from an EoS include the fugacity coefficients/fugacity, residual Helmholtz/Gibbs free energy, residual enthalpy, residual entropy, residual heat capacity, volume/density, and/or their derivatives with respect to the number of moles, the temperature and pressure, etc. The complete relationships between these properties and the residual Helmholtz free energy and its derivatives are described in detail in the book written by Michelsen and Mollerup.22 The most important conceptual property discussed in this work, the reduced residual Helmholtz free energy (F), is defined in eq 1, which will be denoted hereafter as the Helmholtz function: F ¼
Ar ðT, V , nÞ RT
ð1Þ
Generally, the fugacity coefficients and/or their derivatives, with respect to the number of moles, temperature, and pressure are the most important properties when phase equilibria are considered,2325 and can be calculated by the following equation: DF ln Z ln ji ¼ ð2Þ Dni T , V , nj To calculate these properties, the first step is to calculate the corresponding volume for a given temperature, pressure, and number of moles. For general noncubic EoS, it is usually necessary to calculate the first and/or second derivatives of pressure, with respect to total volume, which requires the second and/or third derivatives of the Helmholtz function, with respect to the total volume, depending on the volume/density rootfinding algorithm. Second, with the total volume known, the derivatives of the Helmholtz function with respect to n and then the compressibility factor will be calculated sequentially. The calculation procedure is presented schematically in Figure 1. The calculations of other properties have similar scenarios, except using different combinations of the Helmholtz function and its derivatives. Thus, the entire set of Helmholtz functions, together with their first and second derivatives, with respect to the number of moles, temperature, and total volume, is required in a typical engineering calculation. A class is abstracted to provide the interface to calculate the object function and its derivatives, where each operator is allowed to be called individually, as the class “General Operators” illustrated in Figure 2. The third derivative of the Helmholtz function, with respect to total volume, is also added to this interface, for general purposes. The different algorithms need to be implemented by deriving from this interface class. Another class, called “General Properties”, which provides the access to the object function and its
Figure 1. Schematic procedure to calculate the fugacity coefficients.
derivatives independently, is introduced to match with this General Operators structure. In addition to the applications through its derivatives, the total Helmholtz function can also be approached as a linear combination of separate terms—basic repulsive and attractive contributions plus specific contributions from specific phenomena, such as association, polarity, and electrostatics, etc., as done for the CPA and SAFT familiy EoS.8,9,13,14,17 The total Helmholtz function of CEOS, CPA, and SAFT EoS can be expressed as eq 3, where the CEOS repulsive and attractive contributions are given by explicit formulas,22 and other contributions are implicated by superscript notations. The hard-sphere term is the repulsive contribution, whereas the hard-chain and dispersion terms are the attractive contributions for SAFT. The association term is the specific contribution for both CPA and SAFT.
F total ¼ F rep þ F att þ F SC ( CPA F CEOS þ F Assoc ¼ SAFT F HardS þ F Dispe þ F Chain þ F Assoc
ð3Þ
where B, D, δ1, and δ2 are parameters,22 and “SC” refers to specific contributions. The fact that other properties can be derived from the total Helmholtz function using general relationships and the total Helmholtz function itself is separable provides the foundation for the approach to construct a systematic OO program structure with easy extensibility, reusability, and maintainability. In the OO program architecture, the terminology “HAS-A” refers to a relationship where one object is a part of another object, which behaves according to the rules of ownership.26,27 As the strategy illustrated in Figure 2, the possible Helmholtz function contributions with different physical meanings are implemented separately in a code pool by deriving from General Operators. They can be easily aggregated to a complete thermodynamic package, 6405
dx.doi.org/10.1021/ie101547t |Ind. Eng. Chem. Res. 2011, 50, 6404–6412
Industrial & Engineering Chemistry Research
ARTICLE
Figure 2. Matrix relations between the EOS Package, Specific Packages, XXX_rrA_Opers, and XXX_rrA_Props. [Legend: XXX represents arbitrary symbols; solid boxes with bold text represent class; boxes outlined by dotted lines are class members; boxes with bold and italic text represent explanations; the box outlined by a heavy line represents the Helmholtz Code Pool; boxes outlined by dashed lines are arbitrary classes; lines with arrows represent “derived” relationships; and lines with diamonds represent “composition” relationships.]
called the “Specific Package”, which contains information regarding the number of Helmholtz function terms that it has, what they are, and the necessary routines to connect them. In order to offer a common interface for calls from other components, such as unit operations or thermodynamic analysis utilities in the same process simulator an abstract class, “EOS Package”, could be generalized for all of the Specific Packages. Through this interface class, the properties involved in phase equilibria, as well as thermophysical and volumetric calculations, can be determined, regardless of which Specific Package is used. In the OO program architecture, this relationship can be represented by the terminology “IS-A”,26,27 which means that Specific Package is a specialization of the EOS Package, and EOS Package is a generalization of Specific Package. The general functions of calculating the derived properties from the total Helmholtz function and its derivatives are implemented at the EOS Package level. The pseudocode of calculating the derivative of the total Helmholtz function, with respect to n, is illustrated in the top left of Figure 2, as an example. In addition, the general total volume/density root-finding algorithm, such as a second-order NewtonRaphson method, is also implemented at the EOS Package level. It can be overwritten by the Specific Packages, such as an analytic volume root-finding algorithm for CEOS. Non-CEOS could also have their special procedures for guessing the initial volume and other special treatments, such as how to handle the situation when no corresponding volume root exists for the specified phase status.2830 The pure-compound-related parameters, as well as binary interaction parameters or their coefficients, are read or updated in the Initialize routine. Mixture parameters are calculated in the
Mixingrule routine, where temperature and/or volume dependencies are also considered. In general, Specific Packages will have their own implementations for these two routines. The detailed matrix relationship map of EOS Package, Specific Package, General Operators, and General Properties is illustrated in Figure 2.
3. SPECIFIC PACKAGE STRUCTURE DESIGN: THE PCSAFT EXAMPLE In this section, the PC-SAFT equation of state1315 is taken to illustrate the Specific Package structure design, based on the strategy above. 3.1. PC-SAFT Equation of State and Package Structure. PC-SAFT is a state-of-the-art thermodynamic model with a sound physical basis, which has been applied successfully for an entire variety of substances ranging from small gas molecules up to organic solvents and polymeric systems,1315 and systems with complex intermolecular interactions such as polar interactions, electrical interaction, and induced interactions that only occur in mixtures but not in the pure compounds.3133 The original PC-SAFT version with an association term will be discussed here.13,14 The Helmholtz function of the PC-SAFT EOS can be expressed as
FðT, V , nÞ ¼ ðF HardS þ F Chain Þ þ F Dispe þ F Assoc
ð4Þ
where the superscripts refer to the hard-sphere, hard-chain, dispersion, and association contributions of the total Helmholtz 6406
dx.doi.org/10.1021/ie101547t |Ind. Eng. Chem. Res. 2011, 50, 6404–6412
Industrial & Engineering Chemistry Research
ARTICLE
Figure 3. PC-SAFT Package structure and PC-SAFT Parameters structure. [Legend: boxes outlined by solid lines with bold text represent class; boxes outlined by dotted lines represent class members; the box outlined by solid lines with italic text represents a user-defined simple type (enumeration); lines with arrows represent “derived” relationships; and lines with diamonds represent “composition” relationships.]
function, respectively. Their detailed expressions can be found in the Appendix. As discussed above, they are implemented individually in the code pool. The volume-related parameters ξl, which is defined by eq A.8 in the Appendix, and the radial distribution function (RDF) defined in eq A.11 in the Appendix, are also widely used in other EoS, so it is useful to have general implementations for them as the Helmholtz functions do. This is the reason that the general interface class is called “General Operators” ratehr than “EOS Operators” or “Helmholtz Operators”. The scheme of the PCSAFT Package is illustrated in Figure 3, where the detailed EOS Package members have not been included. It is useful to design a general structure for the PC-SAFT parameters, which can also be used in other EoS, such as CPA and the group-contribution (GC) PC-SAFT.17 To keep its generality, one compound is given the flexibility to have more than one segment. As illustrated in Figure 3, the PC-SAFT parameter class holds an integer member for the number of segment types that it has and a vector member for these segment properties. The segment properties class holds the pure PCSAFT EOS segment parameters—segment number (m), segment diameter (σ), segment energy (ε), and association properties—if it is an associating segment. The association properties class holds an integer member for the number of association types that it has, an integer member for the number of equal sites that the association type has, and an vector member of association site properties, which holds association site type, association
energy, and association volume. The association site type is defined by the atom types, virtual types, or just acceptordonor types, according to engineering preference, if the cross association rule can be defined in advance. From the expressions listed in the Appendix, it could be seen that the hard-sphere, hard-chain, and dispersion terms become based on segment fraction instead of composition, according to this PC-SAFT parameters structure, and the association term still based on site fraction. Therefore, the original PC-SAFT can be considered as a special case of this structure. It is flexible; thus, no code changes are required when the segment types and/or associating schemes8 change. 3.2. Property Calculation Logic: Fugacity Coefficient Example. Taking fugacity coefficient calculation as an example, Figure 4 illustrates the main calculation flowchart of the PCSAFT Package. Updates of the parameters for the Helmholtz functions will be executed at their own level, while the general relationships between derived properties and the total Helmholtz function will be executed at the EOS Package level (see blocks 1 and 2 in Figure 4). Generally, the updates of parameters can be classified into three categories, according to their dependencies: only compound-dependent, temperature-dependent also, and/ or volume-dependent also. They would be further classified for pure compounds and mixture, which are implemented in the Initialize and Mixingrule routines, respectively, as discussed above. These classifications are useful for further code efficiency optimization in large process simulations. 6407
dx.doi.org/10.1021/ie101547t |Ind. Eng. Chem. Res. 2011, 50, 6404–6412
Industrial & Engineering Chemistry Research 3.3. An Example for Package Validation. The binary mixture of water and 1-pentanol is used to illustrate whether the strategy
Figure 4. Fugacity coefficient calculation flowchart for the PC-SAFT Package. Block 1 is executed in the PC-SAFT Package level; block 2 is executed in the EOS Package level. [Legend: Y denotes yes, and N denotes no; solid lines with arrows represent the “flow of control”, and dotted lines with arrows represent “explanations”.]
ARTICLE
works or not. As shown in Figure 5, the phase equilibria behavior reported by others14 is reproduced here. This is because our strategy just rearranges the calculating routines of the EoS, but has nothing to do with their calculation result accuracies. The contributions of the separate terms of the total Helmholtz function can be easily tracked. As shown in Table 1, for this mixture, at vaporliquid equilibrium conditions of T = 383.15 K and P = 101.325 kPa, the separate Helmholtz functions and their derivatives, with respect to the number of moles (n) and total volume are obtained at the EOS Package level, regardless of which Specific Package is used.
4. RESULTS This OO strategy is suitable to build large and growing thermodynamic calculation programs, where many similar parts are used in different EoS. For example, the same association, polar, and electrolytic contributions can be used in different advanced EoS. The schematic hierarchy structure of the entire EoS program is illustrated in Figure 6. It is easy to add a new variant and construct a new package from the existing components; and it is easy for new developers to quickly begin working, and maintain quick continuous developments for the entire program with negligible code efficiency loss. 4.1. Easy Extensibility. If a specific cubic equation of state, such as PR, has been implemented, it can be easily generalized and extended for SRK, PRSV, TST, etc. The general formulas of CEOS can be found elsewhere.22 The different parts of the general parameters and binary interaction parameters are implemented in the Initialize routine. In addition, if the asymmetric binary interaction parameters are introduced to the PRSV EoS, the Mixingrule routine must be implemented in the PRSV package. The general total volume root-finding algorithm is overwritten by the analytical algorithm at the Cubic EOS Package level, as the second layer illustrated in Figure 6. The same logic is adopted for the SAFT family EoS developments, as the third layer illustrated in Figure 6. Only ∼1%5% and 10%50% of the amount of code must be implemented for new variants of the CEOS and SAFT family EoS, respectively. Meanwhile,
Figure 5. Isobaric, heteroazeotropic vaporliquid, and liquidliquid equilibria of a water1-pentanol mixture at a pressure of P = 101.325 kPa. Comparison of experimental data34,35 to calculation results of PC-SAFT (kij = 0.016). 6408
dx.doi.org/10.1021/ie101547t |Ind. Eng. Chem. Res. 2011, 50, 6404–6412
Industrial & Engineering Chemistry Research
ARTICLE
Table 1. Separate Reduced Residual Helmholtz Free-Energy Contributions of a Mixture of Water and 1-Pentanol at T = 383.15 K and P = 101.325 kPa Derivatives wrt Mole Number n Reduced Residual Helmholtz Function (mol)
Derivatives wrt Total Volume (mol/m3)
wrt Compound Water
wrt Compound 1-Pentanol
type
vapor phase
liquid phase
vapor phase
liquid phase
vapor phase
liquid phase
vapor phase
liquid phase
total HardS Chain Dispe Assoc
0.01741 4.7087 103 1.3863 103 0.01191 8.8256 103
5.3049 8.6837 2.5679 10.1167 1.3040
5.6110 104 1.5253 104 4.4885 105 3.8533 104 2.8341 104
10.0792 159.6980 33.9599 119.9459 15.8715
0.03122 4.8236 103 7.5761 104 0.01445 0.02084
4.5408 6.6727 0.9144 7.4250 2.8741
0.04277 0.01989 7.3609 103 0.04514 0.01016
6.7479 28.9890 7.1964 25.6668 2.8737
Figure 6. Schematic structure of the entire EoS program. Boxes outlined by solid lines with bold text represent class; boxes outlined by dotted lines represent class members; and lines with arrows represent “derived” relationships.
Figure 7. Schematic diagram of CPA and SAFT Helmholtz functions composition. Solid lines represent the separate Helmholtz function contribution; dotted lines represent the composition/construction relationship.
developers do not need to take the risk of missing functions, errors from variables renaming, and revalidating routines when transferring code. 4.2. Easy Reusability. If the separate general algorithms of CEOS and the association term of SAFT family EOS already exist in the code pool, it is easy to construct a package for CPA, as sketched in Figure 7. In fact, it is true that the SRK EoS and the original SAFT were developed and applied in industries before CPA. The only work developers must do for CPA Package development is to update the related parameters, for which the PC-SAFT parameters structure can be reused, because they have the same number of pure compound parameters. 4.3. Easy Maintainability. The OO hierarchy structure opens a clear perspective, which allows new developers to get familiar
with the structure and begin their work quickly. It is especially beneficial to the long-term developments for large systems. If each model has its own specific implementations for all parts in a large system, developers may get lost when deep bugs happen. The strategy of putting the general relationships at the EOS Package level allow developers to easily track which parts are not coded correctly, as shown in Table 1, where not all of the derivatives are included, because other properties can be obtained in the same way. 4.4. Program Performance. Program performance is very important in process simulations, especially for dynamic simulations, such as Operator Training Simulator (OTS) and online applications. The significance of our strategy will become very minor if the extensibility, reusability and maintainability are purchased at the expense of code efficiency loss. Ten million times fugacity coefficients calculations for a 10-compound system was used to test the code efficiency for CEOS and PC-SAFT developed by this approach. The comparisons with the mature commercial code are shown in Table 2. It could be seen that the code efficiency loss is
