## Estimating Binary Interaction Parameter for Peng Robinson EOS

This article shows how to estimate binary interaction parameter (BIP) used in Peng Robinson (PR) Equation of State (EOS) from experimental data by regression in excel spreadsheet.

Example
Determine binary interaction parameter used in Peng Robinson EOS for a mixture of acetone and chloroform from Tx-y experimental data available at 101.33 kPa.

Obtain pure component properties like critical temperature (Tc), critical pressure (Pc) and accentric factor ω from literature. Experimental K1 is determined from x, y data.

``K1Experimental = y1 / x1``

Peng Robinson EOS parameters are calculated.

````κi = 0.37464 + 1.54226ω - 0.26992ω²`
`αi = [ 1 + κi (1 - (T/Tc)0.5)]²`
`ai = 0.45724 (RTc)²α / Pc`
`bi = 0.07780 RTc / Pc````

Mixture parameters are calculated next

````a12= [(a1.a2)0.5(1 - k12)] = a21`
`a = a1.x1² + 2.a12.x1.x2 + a2.x2²`
`b = b1.x1 + b2.x2`
`A = aP/(RT)²`
`B = bP/RT````

where k12 is Binary Interaction Parameter (BIP) to be estimated. For first run, use a guess value to start calculation.

Featured Resources :

Hydrocarbon Processing the premier magazine providing job-help information to technical and management personnel in petroleum refining, gas processing, petrochemical/chemical and engineer/constructor companies throughout the world – since 1922.

Following cubic equation is solved to get ZL.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Above equation can be written as

``Z³ + C2.Z² + C1.Z + C0 = 0``

### Solving Cubic Equation

Cubic equation is solved using following procedure. Calculate Q1, P1 & D.

````Q1 = C2.C1/6 - C0/2 - C2³/27`
`P1 = C2²/9 - C1/3`
`D = Q1² - P1³````

If D >= 0, then equation has only one real root provided by

``Z1 = (Q1 + D0.5)1/3 + (Q1 - D0.5)1/3 - C2/3``

If D < 0, then equation has 3 real roots, following parameters are calculated

````t1 = Q1² / P1³`
`t2 = (1 - t1)0.5 / t10.5. Q1/abs(Q1)`
`θ = atan(t2)````

Roots are calculated as following –

````Z0 = 2.P10.5.cos(θ/3) - C2/3`
`Z1 = 2.P10.5.cos((θ + 2*Π)/3) - C2/3`
`Z2 = 2.P10.5.cos((θ + 4*Π)/3) - C2/3````

Roots thus calculated are arranged in descending order, highest root gives ZV and lowest root gives ZL.

### Fugacity

Based on ZL, liquid fugacity φiL is calculated for each data set.

In next step, Vapor phase fugacity is calculated. Mixture properties are estimated as following –

````a = a1.y1² + 2.a12.y1.y2 + a2.y2²`
`b = b1.y1 + b2.y2`
`A = aP/(RT)²`
`B = bP/RT````

Cubic equation is solved as per method shown above to get ZV.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Based on ZV, vapor fugacity φiV is calculated for each data set.

Featured Resources :

Serves chemical engineering professionals in the chemical process industry including manufacturing, engineering, government, academia, financial institutions and others allied to the field serving the global chemical process industry.

K1 is calculated as:

``K1Calculated = φ1L/φ1V``

Difference of K1Experimental and K1Calculated is obtained for all data points. An objective function is defined as summation of all these differences. Click on Solver in Data Ribbon (Excel 2010) to open dialog box for Solver parameters and input data as shown below.

Minimize the objective function by changing values of kij. Uncheck Make Unconstrained Variables Non-Negative, as k12 can take negative value. Click solve to start regression and new value k12 is calculated.

`` k12 = -0.0605``

Above value is then used to calculate y1 values and results are plotted to check the deviation.

Example
Determine temperature dependent binary interaction parameter used in Peng Robinson EOS ( k12 = k1 + k2.T + k3/T ) for a mixture of acetone and n-Hexane from P x-y experimental data available at 318.15 °K.

Follow above steps to calculate pure component and mixture parameters. Change formula for kij to make it temperature dependent. As an initial guess put values for 3 parameters (k1,k2,k3) which gives a numeric difference value. In solver select these 3 variables to be varied for minimizing the difference value and then solve it. After doing regression binary interaction parameters are obtained and result is plotted as following.

````k1 = -0.0218`
`k2 = 3.142x10-4`
`k3 = 10.155````

Note : In case of non-convergence, try different initial guess values for kij and then let solver obtain optimized values.

Spreadsheet for Regressing BIP for Peng Robinson EOS

## PT Flash Calculation using PR EOS

PT Flash calculation determines split of feed mixture F with a molar composition Zi, into Vapor V and Liquid L at pressure P and temperature T. These calculations can be done in a excel spreadsheet using Peng Robinson Equation of State (PR EOS). To start with bubble point pressure (PBubble) and dew point pressure (PDew) are determined for feed mixture.

• P < PDew, Mixture exists as super-heated vapor.
• P > PBubble, Mixture exists as sub-cooled liquid.
• PDew < P < PBubble, mixture exist in vapor and liquid phase.

Initial guess of vapor fraction V and Ki is made as following.

````V = (PBubble - P)/(PBubble - PDew)`
`Ki = exp[ ln(Pc/P) + ln(10)(7/3)(1 + ω )(1-T/Tc)]````

Based on initial Ki values, iteration is done to get value of V which satisfies material balance on system.

````Yi = Ki.Xi`
`1 = V + L`
`Zi = V.Yi + L.Xi````

where V & L are vapor and liquid fractions. Solving above equations for Xi gives :

``Xi = Zi / ( V.( Ki - 1) + 1 )``

Featured Resources :

Hydrocarbon Processing the premier magazine providing job-help information to technical and management personnel in petroleum refining, gas processing, petrochemical/chemical and engineer/constructor companies throughout the world – since 1922.

At Flash conditions

``Σ Yi - Σ Xi = 0``

Above equation can be solved by iteration using Newton Raphson method. Function F(V) is defined as:

````F(V) = Σ Yi  - Σ Xi`
`F(V) = Σ [Zi (Ki - 1)/( V.(Ki - 1) + 1)]````

Derivative of F(V) is calculated as:

``F'(V) = Σ -[Zi(Ki - 1)² /( V.(Ki - 1) + 1)²]``

New estimate of vapor fraction is calculated as:

``V New = V - F(V)/F'(V)``

Function F(V) and F'(V) are calculated based on new vapor fraction and this process is repeated till there is negligible difference in between V and VNew. Vapor fraction thus obtained is then used to estimate vapor and liquid molar composition (Yi & Xi).

### Iteration for Ki

Vapor (Yi) and Liquid (Xi) mol fractions estimated above are used to generate values for Ki. Parameters for Peng Robinson EOS are calculated for each component i.

````κi = 0.37464 + 1.54226ω - 0.26992ω²`
`αi = [ 1 + κi (1 - (T/Tc)0.5)]²`
`ai = 0.45724 (RTc)²α / Pc`
`bi = 0.07780 RTc / Pc````

#### φiL Calculation

Mixture parameters are calculated.

````aij = [(ai.aj)0.5(1 - kij)] = aji`
`a = ΣiΣj aij.Xi.Xj`
`b = Σi bi.Xi`
`A = aP/(RT)²`
`B = bP/RT````

where, kij’s are Binary Interaction Parameter available from literature. Following cubic equation is solved to get ZL.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Roots calculated are arranged in descending order, highest root gives ZV and lowest root gives ZL.

Based on ZL, liquid fugacity φiL is calculated for each component.

Featured Resources :

Serves chemical engineering professionals in the chemical process industry including manufacturing, engineering, government, academia, financial institutions and others allied to the field serving the global chemical process industry.

#### φiV Calculation

Mixture parameters are calculated.

````a = ΣiΣj aij.Yi.Yj`
`b = Σi bi.Yi`
`A = aP/(RT)²`
`B = bP/RT````

Cubic equation is solved to get ZV.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Based on ZV, vapor fugacity φiV is calculated for each component.

Ki is calculated as:

``Ki = φiL/φiV``

New values of Ki thus calculated are again used to estimate V and thereafter Xi & Yi. Iteration is repeated till there is no further change in Ki values. Typically, in 10 iterations change in Ki values become negligible.

All above calculations along with iterative procedure for flash calculation have been provided in below spreadsheet.

Spreadsheet for PT Flash calculation using PR EOS

## Dew T Flash using PR EOS

Dew T flash calculation determine dew point temperature (T) and liquid mol fraction (Xi) for a mixture at given pressure (P) and vapor mol fraction (Yi). These calculations can be performed in excel spreadsheet using Peng Robinson Equation of State (PR EOS).

Estimate temperature T and liquid mol fraction (Xi). T can be estimated as following –

````T = Σ Tisat Xi`
`Tisat = Tc/[ 1 - 3.ln(P/Pc)/(ln(10).(7 + 7ω)) ]````

where Pc, Tc and ω are critical constants for a component i. Liquid mol fraction is estimated as following

````Xi = Yi/Ki`
`Ki = exp[ ln(Pc/P) + ln(10)(7/3)(1 + ω )(1-T/Tc)]````

Featured Resources :

Hydrocarbon Processing the premier magazine providing job-help information to technical and management personnel in petroleum refining, gas processing, petrochemical/chemical and engineer/constructor companies throughout the world – since 1922.

First iteration starts with estimated T and Xi. Parameters for Peng Robinson EOS are calculated for each component i.

````κi = 0.37464 + 1.54226ω - 0.26992ω²`
`αi = [ 1 + κi (1 - (T/Tc)0.5)]²`
`ai = 0.45724 (RTc)²α / Pc`
`bi = 0.07780 RTc / Pc````

Mixture parameters are calculated next

````aij = [(ai.aj)0.5(1 - kij)] = aji`
`a = ΣiΣj aij.Yi.Yj`
`b = Σi bi.Yi`
`A = aP/(RT)²`
`B = bP/RT````

where, kij’s are Binary Interaction Parameter available from literature.
Following cubic equation is solved to get ZV.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Above equation can be written as

``Z³ + C2.Z² + C1.Z + C0 = 0``

### Solving Cubic Equation

Cubic equation is solved using following procedure. Calculate Q1, P1 & D.

````Q1 = C2.C1/6 - C0/2 - C2³/27`
`P1 = C2²/9 - C1/3`
`D = Q1² - P1³````

If D >= 0, then equation has only one real root provided by

``Z1 = (Q1 + D0.5)1/3 + (Q1 - D0.5)1/3 - C2/3``

If D < 0, then equation has 3 real roots, following parameters are calculated

````t1 = Q1² / P1³`
`t2 = (1 - t1)0.5 / t10.5. Q1/abs(Q1)`
`θ = atan(t2)````

Roots are calculated as following –

````Z0 = 2.P10.5.cos(θ/3) - C2/3`
`Z1 = 2.P10.5.cos((θ + 2*Π)/3) - C2/3`
`Z2 = 2.P10.5.cos((θ + 4*Π)/3) - C2/3````

Roots thus calculated are arranged in descending order, highest root gives ZV and lowest root gives ZL.

Featured Resources :

Serves chemical engineering professionals in the chemical process industry including manufacturing, engineering, government, academia, financial institutions and others allied to the field serving the global chemical process industry.

### Fugacity

Based on ZV, vapor fugacity φiV is calculated for each component.

As a next step, Liquid phase fugacity is calculated. Mixture properties are estimated as following –

````a = ΣiΣj aij.Xi.Xj`
`b = Σi bi.Xi`
`A = aP/(RT)²`
`B = bP/RT````

Cubic equation is solved as per method shown above to get ZL.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Based on ZL, liquid fugacity φiL is calculated for each component.

Liquid phase mol fraction is calculated as

``Xi = Yi.φiV/φiL``

New values of Xi thus calculated are again used to estimate φiL and thereafter Xi. This iteration is repeated till there is no further change in Xi values. Typically, in 25 iterations change in Xi values become negligible.

At the end of iteration ΣXi is calculated, if it is close to 1, results are obtained. If not, new value of T is estimated such that ΣXi is close to 1. In excel it can be achieved by using GOAL SEEK function, in which T value is changed to make summation equal to 1.

### Note

For some initial values of Temperature, Xi become equal to Yi and summation ΣXi becomes 1, it happens when initial guess for T falls in critical region. For such cases use different value of temperature, such that summation is not equal to 1 and then use Excel GOAL SEEK function to estimate Dew Point Temperature and liquid mol fractions Xi.

All above calculations along with iterative procedure for flash calculation have been modeled in below spreadsheet.

Spreadsheet for Dew T Flash using PR EOS

## Bubble T Flash using PR EOS

Bubble T flash calculation determine bubble point temperature (T) and vapor mol fraction (Yi) for a mixture at given pressure (P) and liquid mol fraction (Xi). These calculations can be performed in excel spreadsheet using Peng Robinson Equation of State (PR EOS).

Estimate temperature T and vapor mol fraction (Yi). T can be estimated as following –

````T = Σ Tisat Xi`
`Tisat = Tc/[ 1 - 3.ln(P/Pc)/(ln(10).(7 + 7ω)) ]````

where Pc, Tc and ω are critical constants for a component i. Vapor mol fraction is estimated as following

````Yi = Ki Xi`
`Ki = exp[ ln(Pc/P) + ln(10)(7/3)(1 + ω )(1-T/Tc)]````

Featured Resources :

Hydrocarbon Processing the premier magazine providing job-help information to technical and management personnel in petroleum refining, gas processing, petrochemical/chemical and engineer/constructor companies throughout the world – since 1922.

First iteration starts with estimated T and Yi. Parameters for Peng Robinson EOS are calculated for each component i.

````κi = 0.37464 + 1.54226ω - 0.26992ω²`
`αi = [ 1 + κi (1 - (T/Tc)0.5)]²`
`ai = 0.45724 (RTc)²α / Pc`
`bi = 0.07780 RTc / Pc````

Mixture parameters are calculated next

````aij = [(ai.aj)0.5(1 - kij)] = aji`
`a = ΣiΣj aij.Xi.Xj`
`b = Σi bi.Xi`
`A = aP/(RT)²`
`B = bP/RT````

where, kij’s are Binary Interaction Parameter available from literature.
Following cubic equation is solved to get ZL.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Above equation can be written as

``Z³ + C2.Z² + C1.Z + C0 = 0``

### Solving Cubic Equation

Cubic equation is solved using following procedure. Calculate Q1, P1 & D.

````Q1 = C2.C1/6 - C0/2 - C2³/27`
`P1 = C2²/9 - C1/3`
`D = Q1² - P1³````

If D >= 0, then equation has only one real root provided by

``Z1 = (Q1 + D0.5)1/3 + (Q1 - D0.5)1/3 - C2/3``

If D < 0, then equation has 3 real roots, following parameters are calculated

````t1 = Q1² / P1³`
`t2 = (1 - t1)0.5 / t10.5. Q1/abs(Q1)`
`θ = atan(t2)````

Roots are calculated as following –

````Z0 = 2.P10.5.cos(θ/3) - C2/3`
`Z1 = 2.P10.5.cos((θ + 2*Π)/3) - C2/3`
`Z2 = 2.P10.5.cos((θ + 4*Π)/3) - C2/3````

Roots thus calculated are arranged in descending order, highest root gives ZV and lowest root gives ZL.

Featured Resources :

Serves chemical engineering professionals in the chemical process industry including manufacturing, engineering, government, academia, financial institutions and others allied to the field serving the global chemical process industry.

### Fugacity

Based on ZL, liquid fugacity φiL is calculated for each component.

As a next step, Vapor phase fugacity is calculated. Mixture properties are estimated as following –

````a = ΣiΣj aij.Yi.Yj`
`b = Σi bi.Yi`
`A = aP/(RT)²`
`B = bP/RT````

Cubic equation is solved as per method shown above to get ZV.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Based on ZV, vapor fugacity φiV is calculated for each component.

Vapor phase mol fraction is calculated as

``Yi = Xi.φiL/φiV``

New values of Yi thus calculated are again used to estimate φiV and thereafter Yi. This iteration is repeated till there is no further change in Yi values. Typically, in 25 iterations change in Yi values become negligible.

At the end of iteration ΣYi is calculated, if it is close to 1, results are obtained. If not, new value of T is estimated such that ΣYi is close to 1. In excel it can be achieved by using GOAL SEEK function, in which T value is changed to make summation equal to 1.

### Note

For some initial values of Temperature, Yi become equal to Xi and summation ΣYi becomes 1, it happens when initial guess for T falls in critical region. For such cases use different value of temperature, such that summation is not equal to 1 and then use Excel GOAL SEEK function to estimate Bubble Point Temperature and vapor mol fractions Yi.

All above calculations along with iterative procedure for flash calculation have been modeled in below spreadsheet.

Spreadsheet for Bubble T Flash using PR EOS

## Specific Heat Ratio of Real Gas

Equation of state is used to derive variety of thermodynamic properties. This article illustrate calculation of specific heat ratio from Peng Robinson Equation of state.

Example
Calculate specific heat ratio ( γ = Cp/Cv ) for methane gas at 11 Bar & 300 °K. Critical constants for Methane are as following

• Critical temperature, Tc : 190.6°K
• Critical Pressure, Pc : 46.002 bar
• Accentric Factor, ω : 0.008

Ideal gas specific heat constants CpIG = A + B.T + C.T² + D.T³ are as following

• A = 4.5980
• B = 0.0125
• C = 2.86 x 10-6
• D = -2.7 x 10-9

where Cp is in cal/mol-K

Peng Robinson equation of state is defined as

`` P = RT / (V - b)  - a / [V(V + b) + b(V - b)]``

where

```` ac = 0.45723553 R²Tc²/Pc`
` b = 0.077796074 RTc/Pc`
` m = 0.37464 + 1.54226ω - 0.26992ω² `
` a = ac[1 + m(1 - (T/Tc)0.5)]²````

Above equation is translated into polynomial form and solved for values of Z using Newton-Raphson method.

```` Z³ - (1 - B)Z² + Z (A - 2B - 3B²) - (AB - B² - B³) = 0`
` Z = PV/RT`
` A = aP/ (RT)²`
` B = bP/ RT````

Featured Resources :

Hydrocarbon Processing the premier magazine providing job-help information to technical and management personnel in petroleum refining, gas processing, petrochemical/chemical and engineer/constructor companies throughout the world – since 1922.

Following partial derivatives are required for calculating thermodynamic properties. First derivative is obtained by differentiation of P with respect to V at constant T.

```` (δP/ δV)T = -RT/(v - b)² + 2a(v + b)/[v(v + b) + b(v - b)]²`
` (δP/ δV)T = -0.00485 bar/(cm3/mol)````

Second derivative is obtained by differentiation of P with respect to T at constant V.

```` (δP/ δT)V = R/(v - b) - a'/[v(v + b) + b(v - b)]`
` (δa/ δT)V = -mac/[(TTc)0.5(1 + m( 1 - (T/Tc)0.5))]`
` (δP/ δT)V = 0.039 bar/K`
` (δT/ δP)V = 25.814 K/bar````

Third derivative is obtained by differentiation of V with respect to T at constant P.

```` (δV/ δT)P = (R/P)[ T(δZ/δT)P + Z]`
` (δZ/ δT)P = Num / Denom`
` Num = (δA/δT)P (B-Z) + (δB/δT)P(6BZ+2Z-3B²-2B+A-Z²)`
` Denom = 3Z&sup2 + 2(B-1)Z + (A-2B-3B²)````

where,

```` (δA/δT)P = (P/(RT)²)(a' - 2a/T)`
` (δB/δT)P = -bP/(RT²)````

### Calculation of Heat Capacities

Ideal gas heat capacity CpIG is calculated at 300 °K from polynomial equation provided above. Specific heat at constant volume for ideal gas, CvIG is calculated using following relation.

`` CvIG = CpIG - R``

Residual heat capacity at constant volume Cv R is calculated from internal energy U R as following.

```` CvR = (δUR/δT)V`
` UR = [(Ta'-a)/b(8)0.5] ln[(Z+B(1+20.5))/(Z+B(1-20.5))]`
` CvR = [Ta"/b(8)0.5] ln[(Z+B(1+20.5))/(Z+B(1-20.5))]````

where,

`` a" = ac m (1 + m)(Tc/T)0.5/ (2TTc)``

Specific heat capacity at constant pressure and volume is calculated using following equation.

```` CpR = CvR + T(δP/δT)V(δV/δT)P - R`
` Cp = CpIG + CpR`
` Cv = CvIG + CvR````

Specific heat ratio is obtained as :

```` γ = Cp / Cv`
` γ = 1.338````

Spreadsheet for Specific heat ratio calculaton from Peng Robinson EOS

## Dew P Flash using PR EOS

Dew P flash calculation determine dew point pressure (P) and liquid mol fraction (Xi) for a mixture at given temperature (T) and vapor mol fraction (Yi). These calculations can be performed in excel spreadsheet using Peng Robinson Equation of State (PR EOS).

Estimate pressure P and liquid mol fraction (Xi). P can be estimated as following –

````P = 1/ Σ Yi/ Pisat`
`Pisat = exp[ ln(Pc) + ln(10)(7/3)(1 + ω )(1-T/Tc)]````

where Pc, Tc and ω are critical constants for a component i. Liquid mol fraction is estimated as following

````Xi = Yi/Ki`
`Ki = exp[ ln(Pc/P) + ln(10)(7/3)(1 + ω )(1-T/Tc)]````

Featured Resources :

Hydrocarbon Processing the premier magazine providing job-help information to technical and management personnel in petroleum refining, gas processing, petrochemical/chemical and engineer/constructor companies throughout the world – since 1922.

First iteration starts with estimated P and Xi. Parameters for Peng Robinson EOS are calculated for each component i.

````κi = 0.37464 + 1.54226ω - 0.26992ω²`
`αi = [ 1 + κi (1 - (T/Tc)0.5)]²`
`ai = 0.45724 (RTc)²α / Pc`
`bi = 0.07780 RTc / Pc````

Mixture parameters are calculated next

````aij = [(ai.aj)0.5(1 - kij)] = aji`
`a = ΣiΣj aij.Yi.Yj`
`b = Σi bi.Yi`
`A = aP/(RT)²`
`B = bP/RT````

where, kij’s are Binary Interaction Parameter available from literature.
Following cubic equation is solved to get ZV.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Above equation can be written as

``Z³ + C2.Z² + C1.Z + C0 = 0``

### Solving Cubic Equation

Cubic equation is solved using following procedure. Calculate Q1, P1 & D.

````Q1 = C2.C1/6 - C0/2 - C2³/27`
`P1 = C2²/9 - C1/3`
`D = Q1² - P1³````

If D >= 0, then equation has only one real root provided by

``Z1 = (Q1 + D0.5)1/3 + (Q1 - D0.5)1/3 - C2/3``

If D < 0, then equation has 3 real roots, following parameters are calculated

````t1 = Q1² / P1³`
`t2 = (1 - t1)0.5 / t10.5. Q1/abs(Q1)`
`θ = atan(t2)````

Roots are calculated as following –

````Z0 = 2.P10.5.cos(θ/3) - C2/3`
`Z1 = 2.P10.5.cos((θ + 2*Π)/3) - C2/3`
`Z2 = 2.P10.5.cos((θ + 4*Π)/3) - C2/3````

Roots thus calculated are arranged in descending order, highest root gives ZV and lowest root gives ZL.

Featured Resources :

Serves chemical engineering professionals in the chemical process industry including manufacturing, engineering, government, academia, financial institutions and others allied to the field serving the global chemical process industry.

### Fugacity

Based on ZV, vapor fugacity φiV is calculated for each component.

As a next step, Liquid phase fugacity is calculated. Mixture properties are estimated as following –

````a = ΣiΣj aij.Xi.Xj`
`b = Σi bi.Xi`
`A = aP/(RT)²`
`B = bP/RT````

Cubic equation is solved as per method shown above to get ZL.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Based on ZL, liquid fugacity φiL is calculated for each component.

Liquid phase mol fraction is calculated as

``Xi = Yi.φiV/φiL``

New values of Xi thus calculated are again used to estimate φiL and thereafter Xi. This iteration is repeated till there is no further change in Xi values. Typically, in 25 iterations change in Xi values become negligible.

At the end of iteration ΣXi is calculated, if it is close to 1, results are obtained. If not, new value of P is estimated such that ΣXi is close to 1. In excel it can be achieved by using GOAL SEEK function, in which P value is changed to make summation equal to 1.

### Note

For some initial values of Pressure, Xi become equal to Yi and summation ΣXi becomes 1, it happens when initial guess for P falls in critical region. For such cases use lower value of pressure such that summation is not equal to 1 and then use Excel GOAL SEEK function to estimate Dew Point Pressure and liquid mol fractions Xi.

All above calculations along with iterative procedure for flash calculation have been modeled in below spreadsheet.

Spreadsheet for Dew P Flash using PR EOS

## Bubble P Flash using PR EOS

Bubble P flash calculation determine bubble point pressure (P) and vapor mol fraction (Yi) for a mixture at given temperature (T) and liquid mol fraction (Xi). These calculations can be performed in excel spreadsheet using Peng Robinson Equation of State (PR EOS).

Estimate pressure P and vapor mol fraction (Yi). P can be estimated as following –

````P = Σ Pisat Xi`
`Pisat = exp[ ln(Pc) + ln(10)(7/3)(1 + ω )(1-T/Tc)]````

where Pc, Tc and ω are critical constants for a component i. Vapor mol fraction is estimated as following

````Yi = Ki Xi`
`Ki = exp[ ln(Pc/P) + ln(10)(7/3)(1 + ω )(1-T/Tc)]````

Featured Resources :

Hydrocarbon Processing the premier magazine providing job-help information to technical and management personnel in petroleum refining, gas processing, petrochemical/chemical and engineer/constructor companies throughout the world – since 1922.

First iteration starts with estimated P and Yi. Parameters for Peng Robinson EOS are calculated for each component i.

````κi = 0.37464 + 1.54226ω - 0.26992ω²`
`αi = [ 1 + κi (1 - (T/Tc)0.5)]²`
`ai = 0.45724 (RTc)²α / Pc`
`bi = 0.07780 RTc / Pc````

Mixture parameters are calculated next

````aij = [(ai.aj)0.5(1 - kij)] = aji`
`a = ΣiΣj aij.Xi.Xj`
`b = Σi bi.Xi`
`A = aP/(RT)²`
`B = bP/RT````

where, kij’s are Binary Interaction Parameter available from literature.
Following cubic equation is solved to get ZL.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Above equation can be written as

``Z³ + C2.Z² + C1.Z + C0 = 0``

### Solving Cubic Equation

Cubic equation is solved using following procedure. Calculate Q1, P1 & D.

````Q1 = C2.C1/6 - C0/2 - C2³/27`
`P1 = C2²/9 - C1/3`
`D = Q1² - P1³````

If D >= 0, then equation has only one real root provided by

``Z1 = (Q1 + D0.5)1/3 + (Q1 - D0.5)1/3 - C2/3``

If D < 0, then equation has 3 real roots, following parameters are calculated

````t1 = Q1² / P1³`
`t2 = (1 - t1)0.5 / t10.5. Q1/abs(Q1)`
`θ = atan(t2)````

Roots are calculated as following –

````Z0 = 2.P10.5.cos(θ/3) - C2/3`
`Z1 = 2.P10.5.cos((θ + 2*Π)/3) - C2/3`
`Z2 = 2.P10.5.cos((θ + 4*Π)/3) - C2/3````

Roots thus calculated are arranged in descending order, highest root gives ZV and lowest root gives ZL.

Featured Resources :

Serves chemical engineering professionals in the chemical process industry including manufacturing, engineering, government, academia, financial institutions and others allied to the field serving the global chemical process industry.

### Fugacity

Based on ZL, liquid fugacity φiL is calculated for each component.

As a next step, Vapor phase fugacity is calculated. Mixture properties are estimated as following –

````a = ΣiΣj aij.Yi.Yj`
`b = Σi bi.Yi`
`A = aP/(RT)²`
`B = bP/RT````

Cubic equation is solved as per method shown above to get ZV.

``Z³ + (B-1)Z² + (A-3B² -2B)Z + (B³+B²-AB) = 0``

Based on ZV, vapor fugacity φiV is calculated for each component.

Vapor phase mol fraction is calculated as

``Yi = Xi.φiL/φiV``

New values of Yi thus calculated are again used to estimate φiV and thereafter Yi. This iteration is repeated till there is no further change in Yi values. Typically, in 25 iterations change in Yi values become negligible.

At the end of iteration ΣYi is calculated, if it is close to 1, results are obtained. If not, new value of P is estimated such that ΣYi is close to 1. In excel it can be achieved by using GOAL SEEK function, in which P value is changed to make summation equal to 1.

### Note

For some initial values of Pressure, Yi become equal to Xi and summation ΣYi becomes 1, it happens when initial guess for P falls in critical region. For such cases use lower value of pressure such that summation is not equal to 1 and then use Excel GOAL SEEK function to estimate Bubble Point Pressure and vapor mol fractions Yi.