A formulation methodology for multicomponent distillation sequences based on stochastic optimization☆

FeiWang ,Yiqing Luo ,*,Xigang Yuan

1 ChemicalEngineering Research Center,SchoolofChemicalEngineering and Technology,Tianjin University,Tianjin 300072,China

2 State Key Laboratory ofChemicalEngineering,Tianjin University,Tianjin 300072,China

1.Introduction

Distillation is the mostwidely used separation technique in petroleum and chemicalprocess industry,playing an importantrole in energy consumption.Process synthesis,which is considered as the main approach ofenergy saving in engineering,has been applied in distillation process optimization decades ago.Generally,the objective ofsynthesis ofdistillation system is to find the bestdistillation con figuration which separates a given multi-componentfeed into individualcomponentsamong a greatnumber offeasible alternatives in terms ofeconomic competence.Due to the importance ofmulti-component distillation system synthesis in industry,a lot ofapproaches have been proposed and reported to solve the problem systematically,including Heuristics[1–3],evolutionary techniques[4,5],superstructure optimization[6–9]and stochastic methods[10–13].

Mathematically,synthesis ofdistillation system can be formulated as a mixed-integer nonlinear programming problem(MINLP)thatsimultaneously optimizes the structuraland parametric variables.However,the application ofmathematicalprogramming approaches to the synthesis of distillation systemis limited by two dif ficulties inevitably.The firstcomes from the combinatorialfeature of the problem.The feasible distillation con figurations increase rapidly as the number of feed compositions increases.The second comes from the non-convexity of the nonlinear formulation that will cause computational challenges.As a result,the mathematicalprogramming approaches are more suitable for solving smallor moderate sized problems.

Based on adaptive random search,stochastic algorithms such as simulated annealing(SA),in which a search can proceed without direct dependence on the features of nonlinear functions,have been proven to be effective in solving the combinatorialoptimization problem.Besides,stochastic algorithms can converge to global optimal solution with maximum probability.The application of SA approach for the synthesis of distillation system was first reported by Floquet et al.[14].But they did not take continuous variable optimization into consideration.An&Yuan[13]developed a robust SA-based method for the synthesis oflarge-scale distillation systems,allowing the use ofnon-sharp splits and heat integration.However,the non-sharp splits were restricted to the case ofone intermediate componentdistributing simultaneously at the top and bottom ofa column in acceptable quantities.

In order to find an optimal separation solution,the research of separation sequences considering non-sharp splits with multiple distribution components is meaningfulin engineering.Thus,the search space for finding the best solution is enlarged.Agrawal et al.[15–18]presented a mathematicalformulation that could be used to generate only the setofbasic con figurations(including non-sharp splits)and described computationalsearch methods thatwere proven to be usefulin the synthesis ofcon figurations with heat duty requirement.

In this paper,by adopting the SA-based method proposed by An&Yuan[13],the separation sequences,in which non-sharp separations of multiple intermediate distribution components are included,are focused on with the minimum totalannualcost as the evaluation function.According to the characteristics of the optimization problem,the previous coding approach used by An&Yuan[13]is no longer applicable.Hence,a new encoding representation for manipulating separation sequences is developed based on the binary sort tree principle.By the new encoding representation,allpossible distillation con figurations for separating N-component feed mixtures into N nearly product streams can be described.To show the use of this method,we start with the case ofnonsplits thatinclude two intermediate distribution components atmost.

2.Problem Formulation

2.1.Problem statement

The problem addressed in this paper can be expressed as:

Given N-component mixture of known conditions:composition,flow rate,temperature and thermodynamics properties,the distillation system for separating the mixture into N products corresponding to the components,with a separation con figuration consisting of N-1 columns that give the lowest totalannualcost should be found.

The assumptions made for the problem are as follows:

(1)Feed mixture cannot form azeotrope.

(2)The feed is assumed to be saturated liquid and the products ofall columns are removed at their bubble point temperature.

(3)The relative volatility ofeach componentis assumed to be constant.

(4)The re flux ratio(R)for any individualcolumn is set to be 1.02–1.5 times as large as the minimum re flux ratio(Rmin)for a column that performs a non-sharp split.The value of R/Rminis fixed at 1.2 for a column that performs a sharp split.

(5)The recoveries ofcorresponding key components vary within a certain range(e.g.≥98%)for a column which performs a non-sharp split.The recoveries of the key components are fixed at a high value(e.g.98%)for a column that performs a sharp split.

(6)Ifa column performs a non-sharp split,it should be assumed that there are atmost two distribution components between the nonadjacent light and heavy key components.

Assumptions(1)–(3)have been widely adopted in distillation system synthesis problems.Assumption(4)is necessary for the lowest cost when each column operates at its optimalconditions dominated by the re flux ratio.Assumption(5)is from the consideration of the importance ofsimple columns(sharp or non-sharp)for industrialapplications.Assumption(6)is to reduce the complexity ofcombination and calculation of the problem.The majority promising flow sheet structures could be included even though above assumptions may restrict the solution from being a “true”optimum.

According to the foregoing problem statement,the optimization problem includes two kinds ofdecision variables:continuous variables ofindividualcolumns and discrete variables of the flow sheet structure.The former indicate the operating pressure(p),the recoveries of the pair of key components(λLK,λHK)and R/Rmin.In discrete variables,heatintegration con figurations and thermalcoupling structures are excluded and only separation sequence is included.Therefore,the optimization problem can be represented in the form:

where C()is the costmodelfor the system,{si}represents the candidate separation sequence and S is a set including all possible separation sequences.

By analyzing the optimization problem,the basic elements can be summarized as follows:

(1)A superstructure which can include allpossible separation structures in theory.

(2)An MINLP modelincluding structuralconstraints,and mass and energy balance equations.

(3)An economy function to evaluate candidate flow sheetstructure.

(4)An optimization algorithm to solve the model.

In the following sub-sections,the aspects willbe discussed in details.

2.2.Superstructure establishment

2.2.1.New coding representation

Consider an N-componentmixture,and the components{C1,C2,…,Cn,Cn+1,…,CN}are ordered based on their volatilities(i.e.C1is the mostvolatile and CNis the least).Thus,an ordered nature numbers 1,2,…,N can be used as the codes to name components in the feed,as shown in Fig.1.

Fig.1.Codes for naming components.

Different from the coding representation used by An&Yuan[13],a new array[i,j]is used,in which “i”represents the light key component and “j”represents the heavy key component,to represent a separate task.Restricted by assumption(6),“i”and “j”can be described via the following if-then statements:if j=i+1,then the code represents a sharp split.If j=i+2,then the code is a non-sharp split with only one intermediate component.If j=i+3,then the code represents a non-sharp split with two intermediate components.For example,for a mixture of N=4,[1,4]represents a separation task C1C2C3|C2C3C4.C1is the light key component,C4is the heavy key component,C2and C3are distribution components.

In order to divide a multi-component feed into individualcomponents,a separation sequence must include multiple separation tasks.Hence,ST(N-1≤ST≤3N-6)is de fined as the number ofseparation tasks.In an initialsequence generated by ST arrays[i,j]randomly,the codes representing sharp splits mustbe included and the codes representing non-sharp splits can be added or removed when the value of ST can be different.Thus,the optimization problem can be divided into finite sub-problems.But it should be mentioned that if a separation task[i,j](j=i+3)exists in an initialsequence{si}ini,the sequence should be updated.Given an initialsequence{si}ini={[1,4],[2,3],[1,2],[3,4]},the product streams of the first column both contain components C2and C3because the first separation task is C1C2C3|C2C3C4.In other words,C2and C3must be separated twice in the simple columns that only have one feed and two products.The solution is to add[i+1,j-1]in an initialsequence.Therefore,the initialsequence{si}ini={[1,4],[2,3],[1,2],[3,4]}can be updated to{si}ini={[1,4],[2,3],[1,2],[3,4],[2,3]}.Theoretically,the value of ST would not be increased and the solution would not affect the integrity of the superstructure.

According to the principle ofcoding approach,there mustbe a one-toone correspondence between separation sequence and separation process structure.However,the separation task order in{si}inidoes not have any speci fic meanings.Take{si}ini={[2,3],[1,2],[3,4]}and{si}ini={[2,3],[3,4],[1,2]}for examples,their separation sequences are respectively C1C2|C3C4–C1|C2–C3|C4and C1C2|C3C4–C3|C4–C1|C2.The two separation sequences are the same in essence while the representations are different.In addition,it must be considered that the codes can be used to denote anyone in the family of separation tasks that have the same pair of key components.For example,an array[2,3]can be used to representthe separation tasks C2|C3,C1C2|C3,C1C2|C3C4or C1C2|C3…CN,etc.A speci fic separation task can be identi fied according to its position in the sequence.For the above purposes,a binary tree is developed in this work.Speci fic contents willbe introduced in the following section.

2.2.2.Encoding and decoding procedure

Binary tree is one of data structures originally used for storing a large collection of data,providing fast access to the collection when one attempts to find,add or remove data.In addition,from a topology graphics point of view,the binary tree structure is similar to a simple distillation column which has one feed and two products,so it can be applied to simulation of distillation process.

By a binary tree,a one-dimensionalunordered initialsequence{si}inican be converted into an ordered separation sequence{si}.A sorting procedure is needed to identify systematically the separation tasks denoted by the elements in the sequence.

To represent a separation sequence,as illustrated in Fig.2,the following binary tree sorting procedure is proposed:

Step 1.Generate an initialarray code series randomly.An initialsequence would be updated when necessary and the rule is listed above.Step 2.A binary tree structure is constructed based on an initialsequence{si}ini.

To start,the firstarray is sorted as a rootnode.The leftand rightsubnodes are searched according to the following rules:

Firstly,for any separation task si([si,1,si,2]),determine the relationship between s1,1and s1,2based on the root node[s1,1,s1,2].Three cases can be identi fied.

Case 1.s1,2=s1,1+1.If si,2≤s1,1,siis sorted as the left sub-node;if si,1≥ s1,2,siis sorted as the right sub-node;if si,1≤s1,1＜ s1,2≤si,2,the crude separation is infeasible,then make an adjustment and empty the binary tree.

Case 2.s1,2=s1,1+2.If si,2＜s1,2,siis sorted as the left sub-node;if si,1＞ s1,1,siis sorted as the right sub-node;if si,1≤s1,1＜ s1,2≤si,2,the crude separation is infeasible,then make an adjustment and empty the binary tree.

Case 3.s1,2=s1,1+3.If si,1≤s1,1,siis sorted as the left sub-node;if si,1＞s1,1and si,2≥s1,2,siis sorted as the right sub-node;if si,1＞s1,1and si,2＜s1,2,siis sorted as the left sub-node when it appears for the first time and siis sorted as the right sub-node when itappears for the second time.

Secondly,the recursive procedure is continued to the end of the sequence.

Step 3.Based on the constructed binary tree by Step 2,a sequence{si}can be informed by a forward tracing procedure(i.e.visits the root node firstly,then the left sub-node and then the right subnode).Another sequence{kj}is generated by applying a backward tracing procedure(i.e.visits its left sub-tree firstly,then the right sub-tree and then the root node).

Step 4.The corresponding separation network can be identi fied via the following rules:for any array si,if si,2＜si-1,2,the feed stream of code simust be the top product of code si-1.If,on the contrary,si,2≥si-1,2,there is no directlink between the two nodes,{kj}is needed to find the connected “father-code”of si.

By such a method,a sequence{si}which is corresponding to a unique separation sequence can be systematically generated.It must be mentioned that the feasibility of an initialsequence{si}inineed not to be inspected separately.A feasibility validation ofan initialsequence is ongoing when a binary tree is being constructed.Meanwhile,an adjustment should take effect and the binary tree established should be empty when an initialsequence is infeasible.Hence,in this paper,the initialsequence that can construct a binary tree completely is always feasible.However,in the literature of An&Yuan[19],the feasibility of a sequence{si}formed by visiting binary tree may be infeasible.

For the demonstration ofFig.2,in the initialsequence{si}ini,the first code[1,4]should be sorted as a root node which can be described as s1,2=s1,1+3.Therefore,according to Case 3,the second code[2,3]should be sorted as a left node for its first time and the third code[3,4]should be sorted as a right node.For the fourth code[1,2],it should be sorted as a left node on the basis of Case 3,and yet code[2,3]has been there.Thus,the left node[2,3]is identi fied as the root node of left subtree and Case 1 should take effect.Then the code[1,2]willbe sorted as a left sub-node of node[2,3].For the fifth code[2,3],it should be sorted as a rightnode for its second time according to Case 3,and yet code[3,4]has been there.Hence,the code[2,3]willbe sorted as a left sub-node of node[3,4].Thus,a binary tree is completely constructed as shown in Fig.2.Then the sorted sequences{si}and{kj}can be obtained by Step 3.Finally,a corresponding separation network can be identi fied by Step 4.In the sequence{si},the first array code[1,4]is identi fied as the split of feed stream,and the separation task is C1C2C3|C2C3C4.For the second code[2,3],its feed stream must be the top product ofprevious code[1,4]because s2,2＜s1,2,and the splitis identi fied as C1C2|C3.A similar situation occurs to the third code[1,2]that is identi fied as the split of the top product of code[2,3].For the fourth code[3,4],since s4,2＞s3,2,a search in{kj}starts to find the code k4=s4,and k5is recorded.Then a search returns to{si}to find the code s1=k5.So the feed stream ofcode[1,4]mustbe the bottom product ofcode s1.The last code[2,3]is identi fied to split the top product ofcode[3,4]based on above rules.By the identification procedure,a complete separation con figuration of illustrative case ofFig.2 can be described as Fig.3.

Fig.2.Separation sequence representation via binary tree.

Fig.3.A separation con figuration ofan illustrative case.

Itis obvious thatthe separation con figuration in Fig.3 can'tbe created as a basic con figuration.Abasic con figuration is identi fied as a con figuration thatuses(N-1)columns to separate N-component feed into N streams each of which is rich in one of the feed components,where each componentis recovered in only one streamfromonly one location.Rakesh Agrawal[16–18]pointed out that non-basic con figurations which recover atleast one componentor mixture in at least two different locations with more than(N-1)column for an N-component feed,should not be included for a computationally ef ficient search space,especially N is greater than five.Because using more columns in the non-basic con figurations was not signi ficant in improving the total heat duty demand but likely increases the capitalcost.Since the view that non-basic con figurations should be excluded in distillation system is adopted,the separation con figuration in Fig.3 is excluded.

2.2.3.Separation sequence transformation

Based on the above identi fication algorithm,a feasible separation sequence is identi fied.Then,itshould be analyzed how a neighboring structure is generated from current structure in order to get allthe possible separation sequences.Considering SA-based algorithmin which neighboring structure generating is the basic operation,a neighboring separation con figuration should be generated from current separation con figuration with the smallest modi fications as far as possible.For above purposes,a rule interchanging the positions ofa pair ofadjacentseparations proposed by Stephanopoulos and Westerberg[4]is adopted in this work.

In this paper,the neighboring structure generating procedures are different according to different sub-problems.

Situation 1.ST=N-1.The positions ofa pair ofadjacentseparations can be directly interchanged because only sharp splits are included in the system.Details are as follows:

Firstly,a code si([si,1,si,2])is chosen randomly.If si,2＜si-1,2,the positions of siand si-1should be interchanged.In this case,the method is consistentwith“upward”rule in evolutionary improvementprocedure.Conversely,if si,2≥ si-1,2,a connected “father-code”of siis needed to be found by identi fication algorithm and the positions of the two nodes should be interchanged.In this case,the method is consistent with“downward”rule in evolutionary improvementprocedure.

Situation 2.ST=3N-6.A neighboring structure can be generated by directly interchanging the positions ofa pair of adjacent codes that represent the splits with two intermediate components.

Situation 3.N-1＜ST＜3N-6.There are two kinds ofmethods to generate a neighboring structure.The first method is to replace a code representing non-sharp split with a new one.The second one is interchanging the positions of a pair of adjacent separations and the rule has been described in Situation 1.

In theory,allthe possible separation sequences can be generated by above approach.

3.Cost Function

In the presentwork,a totalannualcostincludes the annualoperation cost and annualized investment cost.The operating cost is estimated based on the consumption of the cold and heatutilities,and 8500 h ofoperation per year is also assumed.The conventionalshortcut methods of FUG are applied to column designs.The operating pressure of any individual column is related to the overhead and the bottom compositions,the available hottest and coldest utility temperatures,as well as the critical pressure of the components in the mixture.The range of R/RminandλLK,λHKare described in Assumptions(4)and(5),respectively.The relative volatilities are estimated from the bubble and dew point calculation at the speci fied pressure.The method suggested by Guthrie[20]is applied to the calculation of investment cost ofcolumns.The data provided by Douglas[3]is used to estimate the capitalcost of heat exchangers.An overallheat transfer coef ficient is assumed to be 852 W·m-2·K-1for condensers and 568 W·m-2·K-1for reboilers.A capitalcharge factor of0.1 is used to annualize the equipment cost.

4.SA-based Algorithm

SA originates from statisticalmechanics of solid annealing process,i.e.,melting a solid by increasing its temperature firstly and then cooling the solid so that it can crystallize into a minimum free energy state.SA can be identi fied as a stochastic algorithm in which wrong-way movements during the search for the optimum based on an adaptive acceptance or rejection criterion are allowed.

In this paper,the mixed variables optimization problem is solved by applying the improved simplex simulated annealing algorithm developed by An&Yuan[13].The approach is represented by the flowchart in Fig.4.When an initialsolution is generated,the continuous variables are firstly forced to be optimized while the discrete variables are fixed.If the current solution is accepted based on a criterion,the algorithm would proceed by generating a neighboring discrete solution.Conversely,ifthe currentsolution isrejected butthe corresponding con figuration has notbeen rejected,then itis necessary to decide whetherto generate a new discrete solution or to continue the continuous optimum within current discrete solution.Details on the improved simplex simulated annealing algorithm can be found in An&Yuan[13].

5.Case Study

Fig.4.Flowchart of the ISSA algorithm.

An example of a four-component separation problem is studied to demonstrate the proposed approach.In this case,the k-values,enthalpies and associated thermodynamic properties of the component for allthe products are given by the analyticalexpressions of Reid,Prausnitz,and Poling[21].Besides,the minimum temperature differenceΔtminfor heat exchanger is set to 10 K.The feed is saturated liquid at 0.1 MPa.The feed data is given in Table 1 and the data ofavailable utilities is listed in Table 2.

According to the above approach,for the four components in this case,the separation requires at least three sharp separation tasks and at most three non-sharp separation tasks.Thus,the original problem can bedivided into four sub-problems depending on the number ofnon-sharp separation tasks,as listed in Table 3.Each sub-problem willbe solved independently and the best solution among the four sub-problems is the optimum solution of the example.The calculation software for the new coding procedure is programmed using C++language(Microsoft Visual Studio 2010)on an Intel(R)Core(TM)i3-2100/3.10 GHz computer.For practicalconsideration of CPU time,the scales ofcompression and elongation in simplex algorithm vary with speci fic problems.The optimal con figurations for each sub-problem are shown in Fig.5 and the respective optimization results are listed in Tables 4–7.Table 8 lists the economics parametric ofeach sub-problem.

Table 1 Feed data of the example

Table 2 Utility data of the example

Table 3 Four sub-problems of the example

Fig.5.Optimalcon figurations ofeach sub-problem.

Table 4 Optimalresults of Scheme 1 of the example(ST=3)

Table 5 Optimalresults of Scheme 2 of the example(ST=4)

Table 6 Optimalresults of Scheme 3 of the example(ST=5)

To analyze the results thoroughly,separation sequences with nonsharp splits are superior to those with only sharp splits in energy saving and cost saving,as listed in Table 8.The best non-sharp separation con figuration saves as much as 13.9%cost compared with the optimal sharp separation con figuration and 9.7%energy.Furthermore,nonsharp splits with two distribution components between the nonadjacent light and heavy key components perform better than those with one distribution component.The results are qualitatively similar with those of the literature[18].However,in the literature[18],only energy consumption is compared and the con figuration in Scheme 4 performs best.As listed in Table 8,the operating costin Scheme 4 is signi ficantly higher compared with Scheme 3 even though ithas advantage in energy saving.By analyzing Tables 6 and 7,the bubble-point temperature inScheme 4 is higher than thatin Scheme 3.Thus,a higher quality steam is provided so that the operating cost in Scheme 4 is higher based on the given economic parameters.

Table 7 Optimalresults of Scheme 4 of the example(ST=6)

Table 8 Economics parametric comparisons of the example

Ifwe make a restriction that there is only one distribution component between the non-adjacent light and heavy key components,the originalproblem can be divided into three problems.As a result,the most optimalcon figuration is shown in Scheme 2 and the totalannual costis about213041USDwhich is 5.2%higherthan thatofScheme 3.According to the example,the new coding approach enlarges the scope of alternative separation sequences.Meanwhile,the larger the scope is,the better the optimization result is.

6.Conclusions

In the presentwork,a method ofsynthesizing distillation con figurations is proposed.A new coding procedure is proposed in the form of one-dimensionalarray.Besides,with the application of a binary sort tree approach,a robust flow sheet encoding and decoding procedure are developed.The objective is to obtain the optimaldistillation con figuration with N-1 columns for separating mixtures containing N(N≥4)components into N relatively pure products with a minimum totalannualcost.Finally,an improved simulated annealing approach is adopted to solve the MINLP optimization problem and a case is calculated based on FUG method.In the subsequent work,more persuasive examples willbe calculated to prove the practicability and reliability of the method.Currently,the modelonly takes separation sequences into consideration as discrete variables.The further complex con figurations can be considered,such as heat integration and thermalcoupling.

Nomenclature

[1]R.W.Thompson,C.J.King,Systematic synthesis of separation schemes,AIChE J.18(1972)194.

[2]J.D.Seader,A.W.Westerberg,A combined heuristic and evolutionary strategy for synthesis of simple separation sequences,AIChE J.23(1977)951–954.

[3]J.M.Douglas,ConceptualDesign of Chemical Process,McGrawHill,New York,1988.

[4]G.Stephanopoulos,A.W.Westerberg,Studies in process synthesis II:Evolutionary synthesis of optimalprocess flow sheets,Chem.Eng.Sci.43(1976)195–204.

[5]R.Nath,Motard,Evolutionary synthesis of separation processes,AIChE J.27(4)(1981)578–587.

[6]J.E.Hendry,R.R.Hughes,Generating separation process flow sheets,Chem.Eng.Prog.(1972)68–69.

[7]M.J.Andrecovich,A.W.Westerberg,An MILP formulation for heat integrated distillation sequence synthesis,AIChE J.31(1985)1461–1474.

[8]C.A.Floudas,G.E.Paules,A mixed integrated nonlinear programming formulation for the synthesis of heat integrated distillation sequences,Comput.Chem.Eng.12(1988)531.

[9]H.Yeomans,I.E.Grossman,Nonlineardisjunctive programming models forthe synthesis ofheat integrated distillation sequence,Comput.Chem.Eng.23(1999)1135–1151.

[10]Z.Z.Chen,X.G.Yuan,Synthesis of heat integrated non-sharp distillation sequence,Comput.Appl.Chem.14(4)(1997)247–252(in Chinese).

[11]K.F.Wang,Y.Yuan,P.J.Yao,Synthesis and optimization of heat integrated distillation system using an improved genetic algorithm,Comput.Chem.Eng.23(1988)125–136.

[12]E.Marcoulaki,P.Linke,A.Kokossis,Design of separation trains and reactionseparation networks using stochastic optimization methods,Trans.IchemE A Chem.Eng.Res.Des.79(2001)25–32.

[13]W.Z.An,X.G.Yuan,A simulated annealing-based approach to the optimalsynthesis of heat-integrated distillation sequences,Comput.Chem.Eng.33(2009)199–212.

[14]P.Floquet,L.Pibouleau,S.Domenech,Separation sequenecs synthesis:How to use simulated annealing procedure?Comput.Chem.Eng.18(1994)1141–1148.

[15]R.Agrawal,Synthesis of distillation column con figurations for a multicomponent separation,Ind.Eng.Chem.Res.35(4)(1996)1059–1071.

[16]R.Agrawal,Synthesis ofmulticomponent distillation column con figurations,AIChE J.49(2)(2003)379–401.

[17]A.Giridhar,R.Agrawal,Synthesis ofdistillation con figurations:I.Characteristics ofa good search space,Comput.Chem.Eng.34(2010)73–83.

[18]A.Giridhar,R.Agrawal,Synthesis ofdistillation con figurations.II:A search formulation for basic con figurations,Comput.Chem.Eng.34(2010)84–95.

[19]W.Z.An,X.G.Yuan,Synthesis of thermalcoupling complex distillation system based on stochastic optimization algorithm:I.Modeling method,J.Chem.Ind.Eng.57(7)(1997)1591–1598(in Chinese).

[20]K.M.Guthrie,Capitalcost estimating,Chem.Eng.24(1969)114–142.

[21]Reid,R.C.,Prausnitz,J.M.,&Poling,B.E.,The Properties of Gases and Liquids(4th Ed.),New York:McGraw-hill.