CHAPTER 2

CUBIC EQUATION OF STATE

Cubic

equations of state are equations, which when expanded have volume terms raised

to the first, second, and third power. Most commonly encountered phase

equilibrium calculations, such as vapour-liquid equilibria, involve only two

phases for which a cubic equation is suitable. Cubic equations have the

advantage that the three values of volume can be obtained analytically without

the need for an iterative solution procedure. The van der Waals equation of

state (1873) is the simplest cubic equation of state for fluid phase

equilibria. It can be regarded as a “hard sphere term + attractive term”

equation of state composed from the contribution of repulsive and attractive

intermolecular interactions 2. The van der Waals equation was the

first equation capable of representing vapor-liquid coexistence. The pressure (p)

is related to the temperature (T), ideal gas constant (R) and

molar volume (V) via:

……… (2.1)

It has two pure component parameters a and

b. The parameter a is a measure of the attractive forces between

the molecules, and b is related to the size of the molecules. Since the

van der Waals equation of state is cubic in volume, three volumes exist for any

given temperature and pressure. Equation for a and b can be given as:

……… (2.2)

CHAPTER 3

MODELING

There are three basic categories of VLE

calculations which are extensively used in design:

Estimation of equilibrium

temperature, T and vapor or liquid phase compositions, y or x, at constant

pressure (bubble point or dew point temperature calculations).

The prediction of the system

pressure, P, and vapor or liquid compositions, y or x under isothermal

conditions (bubble and dew point pressure calculations).

The evaluation of the phase

compositions, x, y or K-values at a given temperature and pressure (flash

calculations).

Figure 3.1 Basic categories of VLE calculations

In this course project work I used bubble

P calculation to find the Vapor liquid equilibria using ? – ? approach. Steps

to perform Bubble P calculation using three parameter Tsai & Jan EOS to

find vapor phase non-ideality and Wilson equation for liquid phase

non-ideality.

3.1

Bubble

P calculation

During the calculation of the bubble P,

data required are:

Liquid

phase composition at a particular temperature,

Temperature,

Critical

Properties (TC, PC, ZC) of all component that

are present in the system,

Values

of acentric factors of all component that are present in the system,

Antoine

equation constants,

Value

of Wilson parameters for the given system.

Following are the steps for performing

bubble P calculation

1

List

the data required.

2

Find

Vapor pressure of the compounds at the given Temperature using Antoine Equation.

Antoine equation is given as:

……… (3.1)

Where, A, B and C are Antoine constants And Psat

is saturated vapor pressure in mmHg

and T is temperature in Kelvin.

3

Find

activity coefficient by using Wilson equation. Wilson Equation can be given as,

……… (3.2)

……… (3.3)

Where, ? is activity coefficient, A12 and A21 are Wilson equation constants.

4

Use

modified Roult’s to find Bubble P, Modified Roult’s low can be give as,

………

(3.4)

Where, P is total pressure, y is vapor composition, x is liquid composition, ? is fugacity coefficient and ? is activity coefficient and Psat is saturated vapor

pressure.

5

First

assume ? = 1, that reduces equation (3.4) as follow,

……… (3.5)

6

Find

valve of fugacity coefficients by use of Tsai and Jan EOS. The maximum value of fugacity coefficient is

related to vapor phase non-ideality. In this calculation value of P is required that is not available at

start so we have to make an estimate by performing steps 3, 4 and 5.

7

Convert

equation (3.4) in terms of y to find

vapor composition,

……… (3.6)

8

Use

y values found by step 7 in equation (3.4) with ? value that is

found by step 6 to find new value of P.

9

Recalculate

steps 6, 7 and 8 till following

condition is satisfied,

……… (3.7)

10

To

obtain P-x-y data perfume all steps

for different values of x.

3.2

Bubble

T calculation

During the calculation of the bubble T,

data required are:

Liquid

phase composition at a particular temperature,

Pressure,

Critical

Properties (TC, PC, ZC) of all component that

are present in the system,

Values

of acentric factors of all component that are present in the system,

Antoine

equation constants,

Value

of Wilson parameters for the given system.

Following are the steps for performing

bubble T calculation

1

List

the data required.

2

Estimate

value of temperature for the compounds at the given pressure using Antoine

Equation. From equation (3.1) find value of T1sat

and T2sat and

estimated T becomes,

……… (3.8)

3

Find

P1sat and P2sat using Test in equation (3.1),

4

Find

activity coefficient by using Wilson equation. Eq. (3.2) and (3.3),

5

Find

valve of fugacity coefficients by use of Tsai and Jan EOS. The maximum value of fugacity coefficient is

related to vapor phase non-ideality. In this calculation value of T is required that is not available at

start so we have to make an estimate by performing steps 2.

6

Find

Value of y1 and y2 by using Eq. (3.6).

7

Use

y values found by step 6 in equation

(3.4) with ? value that is found by step 5

to find new value of P.

8

Recalculate

steps 6, 7 and 8 till following

condition is satisfied,

……… (3.9)

9

When

above condition satisfied value of temperature obtain is Bubble T. To obtain T-x-y data perfume all steps for different values of x.

3.3

Dew

T calculation

During the calculation of the dew T, data

required are:

Vapor

phase composition at a particular temperature,

Pressure,

Critical

Properties (TC, PC, ZC) of all component that

are present in the system,

Values

of acentric factors of all component that are present in the system,

Antoine

equation constants,

Value

of Wilson parameters for the given system.

Following are the steps for performing Dew

T calculation

1

List

the data required.

2

Estimate

value of temperature for the compounds at the given pressure using Antoine

Equation. From equation (3.1) find value of T1sat

and T2sat and

estimated T becomes,

……… (3.10)

3

Find

P1sat and P2sat using Test in equation (3.1),

4

Take

? = 1 and ? = 1 and use Roult’s law to find x1 and x2,

5

Find

valve of fugacity coefficients by use of Tsai and Jan EOS. (The maximum value of fugacity coefficient is

related to vapor phase non-ideality.)

6

Find

the value of ? by using (3.2) and (3.3).

7

Recalculate

value of P using equation (3.4).

8

Recalculate

steps 2 to 7 till following condition is satisfied,

……… (3.11)

9

When

above condition satisfied value of temperature obtain is Bubble T. To obtain T-x-y data perfume all steps for different values of x.

3.4

Dew

P calculation

During the calculation of the dew P, data

required are:

Vapor

phase composition at a particular pressure,

Temperature,

Critical

Properties (TC, PC, ZC) of all component that

are present in the system,

Values

of acentric factors of all component that are present in the system,

Antoine

equation constants,

Value

of Wilson parameters for the given system.

Following are the steps for performing Dew

P calculation

1

List

the data required.

2

Find

P1sat and P2sat using T in equation (3.1),

3

Take

? = 1 and ? = 1 and use Roult’s law to find x1 and x2,

4

Find

valve of fugacity coefficients by use of Tsai and Jan EOS. (The maximum value of fugacity coefficient is

related to vapor phase non-ideality).

5

Find

the value of ? by using (3.2) and (3.3).

6

Recalculate

value of P using equation (3.4).

7

Recalculate

steps 2 to 7 till following condition is satisfied,

……… (3.12)

8

When

above condition satisfied value of temperature obtain is Bubble T. To obtain T-x-y data perfume all steps for different values of x.

3.5

Flash

calculation

Flash calculation is used when any of the

two values (From P, T, x,

y) are given. In flash calculation K-

value are used for calculation. It does serve as a measure of the “Lightness”

of a constituent species, i.e., its tendency to favor the vapor phase. When its

value is greater that one means species i

exhibits a higher concentration in vapor phase and if less than higher

concentration in liquid phase. K-values can be found from the nomographs.

……… (3.13)

……… (3.14)

……… (3.15)

Above equations are use to for VLE

calculation. Equation (3.13) used when T

& P or x & y values are

given. Equation (3.14) use for bubble point calculation and (3.15) use for dew

point calculation.

Nomographs and Algorithms for all the calculation

are provided in the Appendix.

CHAPTER 4

DESIGN IN EXCEL

In this chapter, the procedure for

performing calculations for bubble P for generating P-x-y data in Excel is

shown. The work in this course project has been restricted to only bubble P

calculation. Following are the steps to perform Bubble P calculation.

Step 1.

Generation

of Proper Data Base System.

Proper Data Base will

make work easy when we need to call values for particular systems. In this work,

I include three different types of the Data Base systems for my calculation.

1.

For

critical properties of individual compounds,

2.

For

binary systems and their Wilson parameters,

3.

For

comparison with experimental data.

First Data Base contains

critical properties, Antoine constants, acentric factor etc. for the

calculation purpose (Fig. 4.1). There are two columns containing no, first for calling the compounds and

second for calling the values.

Figure 4.1 Data

Base for critical properties of components

Second Database (Fig 4.2)contain

Wilson parameter constants and slope and intersects values to generate parity

plot. First column is for giving unique identity for particular system. And

second and third columns are unique identities of first and second compound in

systems respectively. Last two columns are slope and intercepts to generate

parity plot. Ex no is unique identity for the experimental data and use to call

experimental data from Third Data Base. Third Data Base contains P-x-y data of

experiments.

Figure 4.2 Data

Base for Wilson parameter constants

Step 2.

Listing

properties and parameters,

List all the Parameters

and properties that are remain constant during the calculation. For this

calculation following are the parameters and properties that are constant

during the calculation of particular system.

Table 4.1 Constant

Parameters and Properties.

To call the values of different properties

and parameters vlookup function is

use. Syntax for vlookup can be given

as follow,

= vlookup (lookup_value, table_array, column_index_no,

range_lookup)

Step 3.

Design

Tools,

In this work, the button

and combo box from the developer tab were used. To use that one first need to

make visible the developer tab which is hide default in any MS word. To make it

visible follow step below for MS Office 2016 or Office 2013

1.

Go

to File ? Options ? Customize Ribbon,

check the box before the developer.

Figure 4.3 Excel Option for developer tab

2.

From

developer select insert ? combo box,

Figure 4.4 Combo box

3.

After

selecting combo box draw box of appropriate size, and then follow window will

open,

Figure 4.5 Format control for combo box

In input range

select the cells whose value one want to show in box, in cell link select any blank cell, it will show the number of a

system that is selected from the list i.e. if u select fifth system from the

list (combo box) then it will show 5

in the cell and that is helpful to call the values using vlookup.

To create button, go to developer ? insert ? button,

Figure 4.6 Button tool

After selecting button draw a button of an appropriate

size and after that on window will open to ask for macro but we don’t make it

one so just close the window.

Step 4.

Use

the equations mention in chapter 3 for Tsai and Jan to calculate the ? values and Wilson equation to find ?. To calculate follow the step given

bellow,

1.

List

the values of x at uniform intervals,

2.

Find

the values of the parameters required in the calculation of Tsai & Jan EOS

(Calculate till A, B, C).

3.

From

A, B, C generate the equation for Z using equation (2.31). and compare it

with following equation,

………. (4.1)

By doing so,

……… (4.2)

To solve this 3rd

degree polynomial equation following equations are use,

f = ((3*c/a) – (b^2/a^2))/3

g = ((2*b^3/a^3) – (9*b*c/a^2)

+ (27*d/a))/27

h = (g^2/4) + (f^3/27)

i = SQRT ((g^2/4)-h)

j = I ^ (1/3)

k = ACOS(-(g/(2*i)))

L = j*(-1)

M = COS(k/3)

N = SQRT (3) * SIN(k/3)

P = (b/(3*a)) * (-1)

R = -(g/2) + SQRT(h)

S = R ^ (1/3)

T = – (g/2) – SQRT(h)

U = T ^ (1/3)

For h > 0, one real

root and two imaginary

Z1 = (S+U) –

(b/3a),

Z2 = (-

(S+U)/2-(b/3a)) + i (S-U) *SQRT

(3)/2,

Z3 = (-

(S+U)/2-(b/3a)) – i (S-U) *SQRT (3)/2,

For h ? 0, three unique

and real roots

Z1 = 2*j*COS(k/3)

– (b/(3*a)),

Z2 = L*(M+N) +P,

Z3 = L*(M-N) +P;

For f, g, h = 0, only one

real root

Z = (d/a) ^ (1/3) * (-1);

For the all possible

values of Z obtain highest real root

is consider to calculate ? for both

the component from the system and same method is used to find ?b for equation

(2.26).

4.

Use

equation (3.2) and (3.3) to find ?.

5.

Calculate

P1sat and P2sat from Antoine

equation

Step 5.

Bubble

P calculation,

To calculate bubble P in

excel follow the bellow steps,

1.

Make

a column for estimate P, which will

be found by equation (3.5),

2.

Calculate

yi using equation (3.6),

3.

Make

a column for calculated P, which is

calculated by equation (3.4),

4.

Copy

calculated value on estimate as paste as

value only,

5.

Make

column for ?P and ?y. ?P

is difference between estimate and calculated P and ?y is summation of

y values.

6.

Use

solver to make ?y = 1 and ?P = 0,

Step 6.

Make

use of solver,

Solve is not installed by

default in Excel. So, to install and use solver use steps given bellow for

office 2013 or above,

1.

Go

to File ? options ? add-ins ? solver add-in, click on GO.

2.

To

use solver, go to data ? solver,

following window will open,

Figure 4.7 Solver window

3.

Make

set cell as ?P, to value 0 and by

changing cell as copied P cell which

is paste as value only. Add constrain that ?y

cell becomes 1. Run the solver.

Step 7.

Macro

for solver,

Repetitive use of solver

is very tedious and time consuming work. So, to avoid that we have to create macro.

Before that one has to make sure that its macro is enable. For that go to developer ? macro security ? enable all

macro. After that follow steps given bellow,

1.

Go

to developer and click on record macro,

2.

Perform

on single solver process,

3.

Go

to developer and click on stop recording,

4.

Go

to developer and click on Visual Basic,

5.

Copy

the code and paste it in the sheet where calculation is performed,

6.

Sheet

can be found from the list of all active sheet shown in top-left penal,

7.

Go

to tools ? reference and tick mark solver,

press OK,

8.

Edit

the code as follow,

9.

After

creating macro assign it to button mention in step 3.

Sub Macr()

‘

‘ Macr Macro

Range(“BK51:BK71”).Select

Selection.Copy

Range(“BD51″).Select

Selection.PasteSpecial

Paste:=xlPasteValues, Operation:=xlNone, SkipBlanks _

:=False,

Transpose:=False

solverReset

SolverOk

SetCell:=”$BN$51″, MaxMinVal:=3, ValueOf:=0,

ByChange:=”$BD$51″, _

Engine:=1,

EngineDesc:=”GRG Nonlinear”

SolverAdd

CellRef:=”$BO$51″, Relation:=2, FormulaText:=”1″

SolverSolve True

In above coding first

paragraph is add for copy the calculated value and paste as value only so that

solver can be run. In second paragraph first and last line are added fist line

is for reset solver after every titration and last line is to avoid clicking OK every

time solver run. Coding between that is generated by recording macro. Copy

second paragraph and paste it for the number of time one would run it for

different rows or cells. Change the value of cell references with appropriate

values (i.e. if one want to use for next row then $BN$51 becomes $BN$52 and so

on…). Save coding and save excel file as macro enable file.

Step 8.

Graph

plot.

After generating P-x-y data draw a scatter graphs of

1.

P vs x,

2.

P vs y,

3.

y vs x.

Step 9.

Comparison

plots,

To create comparison plot

we have to generate P-x-y data for the

same value of x that is given in experimental data and compare the calculated y and P values with experimental one. After that, draw various comparison

graphs. (All graphs are mentions in appendix.)

CHAPTER 5

RESULTS AND DISCUSSION

In this excel file

around 25 random systems are included and the result obtained are compared with

the actual experimentally 4 generated data. For Ethanol – Benzene

system at 298.15 K temperature following data are generated and comparison with

actual data generate graphs shown in figure 5.1.

Table 5.1 P-x-y

data for Ethanol – Benzene system

x

0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.90

1.00

y

0.00

0.41

0.53

0.59

0.63

0.65

0.67

0.69

0.70

0.74

1.00

P (mmHg)

59.1

90.9

107.5

116.4

121.1

123.6

124.7

125.0

124.4

121.3

95.2

Figure 5.1

Comparison plots for Ethanol – Benzene system

Just like this

more systems can be calculated in Excel. For around 25 azeotropic and non

– azeotropic system calculation has been

done and results are compiled in table 5.2 and table 5.3 with error in pressure,

y values, parity plot error etc.

Figure 5.2 P-x-y diagram and x-y diagram for ethanol –

benzene system

Table 5.2 Result

for non azeotropic systems

Table 5.3 Result

for azeotropic systems

In above tables values of

?Y and ?P average error are found by using equation given bellow,

The number in the final

column are show how good the fitting of the curve is, and 1. Excellent, 2.

Good, 3. Moderate, 4. Poor, 5. Very Poor. Rating is done on the bases of how

many points are deviating from the experimental values and how much is the

deviation. i.e., for benzene – toluene systems all the points are perfectly

match except one so, there is possibility that the point is due to experimental

error, thus it is considered as perfect match.

CONCLUSION

From the above result I

can conclude that,

v

some

of the non azeotropic systems with few of the azeotropic systems shows the very

good match with experimental data, but most of the azeotropic systems are

showing deviation from the experimental values.

v More over system with combination of

aliphatic and aromatic compounds are showing more deviation compare to same

systems.

More details about this

EOS can be obtain by including more compounds.