Flash Calculations using the 
Peng Robinson Equation of State
Authors' Data: Zakia Nasri and Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :
This problem was presented in the Master's report of R. Privat under the direction of J-P Corriou.

Off[General::"spell1"]

Reduced Pressure and Temperature

Tre[i_] := T/Tc[i]

Pre[i_] := P/Pc[i]

Equilibrium constants using the ϕ-ϕ methode and the Peng Robinson Equation of State. Binary interaction parameters for the Peng Robinson EOS were taken from HYSYS 3.2.

ap[i_] := 
      0.45724 a[i]  R^2 Tc[i]^2/Pc[i]

Ap[i_] := ap[i] P/(R T)^2;

bp[i_] := 0.07780 R Tc[i]/Pc[i]

Bp[i_] := P bp[i]/(R T)

(Subscript[ϕ, v])[i_] := 
      Exp[((Subscript[Z, v]) - 1) Bp[i]/(Subscript[B, v]) - Log[(Subscript[Z, v]) - (Subscript[B, v])] - 
          (Subscript[A, v])/(2 (Sqrt[2]) 
                  (Subscript[B, v])) (2 (∑_(j = 1)^4) y[j] Ab[i, j]/
                        (Subscript[A, v]) - 
                Bp[i]/(Subscript[B, v])) Log[((Subscript[Z, v]) + (1 + (Sqrt[2])) 
                      (Subscript[B, v]))/((Subscript[Z, v]) + (1 - (Sqrt[2])) (Subscript[B, v]))]]

(Subscript[ϕ, l])[i_] := 
      Exp[((Subscript[Z, l]) - 1) Bp[i]/(Subscript[B, l]) - Log[(Subscript[Z, l]) - (Subscript[B, l])] - 
          (Subscript[A, l])/(2 (Sqrt[2]) 
                  (Subscript[B, l])) (2 (∑_(j = 1)^4) x[j] Ab[i, j]/
                        (Subscript[A, l]) - 
                Bp[i]/(Subscript[B, l])) Log[((Subscript[Z, l]) + (1 + (Sqrt[2])) 
                      (Subscript[B, l]))/((Subscript[Z, l]) + (1 - (Sqrt[2])) (Subscript[B, l]))]]

a[i_] := (1 + m[i] (1 - Tre[i]^0.5))^2

m[i_] := 
      0.37464 + 1.54226 ω[i] - 0.26992 ω[i]^2

(Subscript[A, v]) = (∑_(i = 1)^4)(∑_(j = 1)^4) y[i] y[j] Ab[i, 
                j];

Ab[i_, j_] = (Ap[i] Ap[j])^0.5 (1 - k[i, j]);

(Subscript[B, v]) = (∑_(i = 1)^4) y[i] Bp[i];

(Subscript[A, l]) = (∑_(i = 1)^4)(∑_(j = 1)^4) x[i] x[j] Ab[i, 
                j];

(Subscript[B, l]) = (∑_(i = 1)^4) x[i] Bp[i];

K[i_] := (Subscript[ϕ, l])[i]/(Subscript[ϕ, v])[i];

k[1, 1] = 0; k[2, 2] = 0; k[3, 3] = 0; 
    k[4, 4] = 0;
k[2, 1] = k[1, 2]; k[3, 1] = k[1, 3]; k[4, 1] = k[1, 4]; 
    k[2, 3] = k[3, 2]; k[3, 4] = k[4, 3]; 
    k[2, 4] = k[4, 2];
k[1, 2] = 0.00185; k[1, 3] = 0.00464; k[1, 4] = 0.00811; 
    k[3, 2] = 0.00132; k[4, 2] = 0.00317; k[4, 3] = 0.00098;

R = 1.987;

Critical Pressure and Temperature and acentric factor for Ethane, Propane, n-Butane and n-Pentane

Pc[1] = 709.8; Tc[1] = 550.0; ω[1] = 0.1064;

Pc[2] = 617.4; Tc[2] = 665.9; ω[2] = 0.1538;

Pc[3] = 550.7; Tc[3] = 765.3; ω[3] = 0.1954;

Pc[4] = 489.5; Tc[4] = 845.6; ω[4] = 0.2387;

Flash Pressure and Temperature are 197.33 °F and 20 atm. Feed is composed of 2.5% Ethane, 35% Propane, 60% n-Butane and 2.5% n-Pentane.

T = 197.33 + 459.67;

P = 20 14.696;

z[1] := 0.025; z[2] = 0.35; z[3] = 0.60; z[4] = 0.025;

Rachford and Rice Equation and Flash Calculations.

Eq[i_] := y[i] ==  K[i] x[i]

V = (Subscript[ϕ, vap]) F; L = (1 - (Subscript[ϕ, vap])) F;

EQ[i_] := 
        x[i] ==  
          z[i]/(1 + (Subscript[ϕ, vap]) (K[i] - 1));

FindRoot[{Eq[1], Eq[2], Eq[3], (∑_(i = 1)^4) y[i] == 1, 
          EQ[1], EQ[2], 
          EQ[3], (∑_(i = 1)^4) x[i] == 
            1, (∑_(i = 1)^4)(( z[i] (K[i] - 1))/(1 + (Subscript[ϕ, vap]) (K[i] - 1))) == 0, 
          (Subscript[Z, l])^3 + (-1 + (Subscript[B, l])) (Subscript[Z, l])^2 + 
              (Subscript[Z, l]) ((Subscript[A, l]) - 3  (Subscript[B, l])^2 - 2  (Subscript[B, l])) - (Subscript[A, l]) (Subscript[B, l]) + (Subscript[B, l])^2 + 
              (Subscript[B, l])^3 == 0, 
          (Subscript[Z, v])^3 + (-1 + (Subscript[B, v])) (Subscript[Z, v])^2 + 
              (Subscript[Z, v]) ((Subscript[A, v]) - 3  (Subscript[B, v])^2 - 2  (Subscript[B, v])) - (Subscript[A, v]) (Subscript[B, v]) + (Subscript[B, v])^2 + 
              (Subscript[B, v])^3 == 0}, {x[1], 0.01, 0.05}, {x[2], 0.3, 
          0.4}, {x[3], 0.6, 0.7}, {x[4], 0.05, 0.1}, {(Subscript[ϕ, vap]), 0.001, 
          0.9}, {y[1], 0.01, 0.05}, {y[2], 0.4, 0.5}, {y[3], 0.5, 
          0.6}, {y[4], 0.01, 0.05}, {(Subscript[Z, l]), 0.07, 0.09}, {(Subscript[Z, v]), 0.65, 0.85}, 
        MaxIterations -> 10000] // Chop

FindRoot[{Eq[1], Eq[2], Eq[3], (∑_(i = 1)^4) y[i] == 1, 
          EQ[1], EQ[2], 
          EQ[3], (∑_(i = 1)^4) x[i] == 
            1, (∑_(i = 1)^4)(( z[i] (K[i] - 1))/(1 + (Subscript[ϕ, vap]) (K[i] - 1))) == 0, 
          (Subscript[Z, l])^3 + (-1 + (Subscript[B, l])) (Subscript[Z, l])^2 + 
              (Subscript[Z, l]) ((Subscript[A, l]) - 3  (Subscript[B, l])^2 - 2  (Subscript[B, l])) - (Subscript[A, l]) (Subscript[B, l]) + (Subscript[B, l])^2 + 
              (Subscript[B, l])^3 == 0, 
          (Subscript[Z, v])^3 + (-1 + (Subscript[B, v])) (Subscript[Z, v])^2 + 
              (Subscript[Z, v]) ((Subscript[A, v]) - 3  (Subscript[B, v])^2 - 2  (Subscript[B, v])) - (Subscript[A, v]) (Subscript[B, v]) + (Subscript[B, v])^2 + 
              (Subscript[B, v])^3 == 0}, {x[1], 0.01}, {x[2], 0.3}, {x[3], 
          0.6}, {x[4], 0.05}, {(Subscript[ϕ, vap]), 0.3}, {y[1], 0.01}, {y[2], 
          0.4}, {y[3], 0.5}, {y[4], 0.01}, {(Subscript[Z, l]), 0.07}, {(Subscript[Z, v]), 0.65}, 
        MaxIterations -> 10000] // Chop

True Vapor Pressure (TVP) Computation using the Peng Robinson Equation of State. We find TVP=108.476 psia.

P =.

(Subscript[ϕ, vap]) = 0;

T = 100 + 459.67;

FindRoot[{Eq[1], Eq[2], Eq[3], (∑_(i = 1)^4) y[i] == 1, 
          EQ[1], EQ[2], 
          EQ[3], (∑_(i = 1)^4) x[i] == 
            1, (∑_(i = 1)^4)(( z[i] (K[i] - 1))/(1 + (Subscript[ϕ, vap]) (K[i] - 1))) == 0, 
          (Subscript[Z, l])^3 + (-1 + (Subscript[B, l])) (Subscript[Z, l])^2 + 
              (Subscript[Z, l]) ((Subscript[A, l]) - 3  (Subscript[B, l])^2 - 2  (Subscript[B, l])) - (Subscript[A, l]) (Subscript[B, l]) + (Subscript[B, l])^2 + 
              (Subscript[B, l])^3 == 0, 
          (Subscript[Z, v])^3 + (-1 + (Subscript[B, v])) (Subscript[Z, v])^2 + 
              (Subscript[Z, v]) ((Subscript[A, v]) - 3  (Subscript[B, v])^2 - 2  (Subscript[B, v])) - (Subscript[A, v]) (Subscript[B, v]) + (Subscript[B, v])^2 + 
              (Subscript[B, v])^3 == 0}, {x[1], 0.01}, {x[2], 0.3}, {x[3], 
          0.6}, {x[4], 0.05}, {P, 100}, {y[1], 0.1}, {y[2], 0.5}, {y[3], 
          0.3}, {y[4], 0.01}, {(Subscript[Z, l]), 0.07}, {(Subscript[Z, v]), 0.85}, 
        MaxIterations -> 1000000] // Chop

Conclusion: Both flash and TVP calculations are in agreement with results of Privat and HYSYS 3.2.