Adiabatic Flash Calculations using the 
Peng Robinson Equation of State
Author's Data: 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 excellent book by Henley and Seader, Equilibrium-Stage Separation Operations in Chemical Engineering, Wiley, 1981.

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.202000007033348; k[1, 3] = 0.285100013017654; 
    k[1, 4] = 0.285100013017654; k[3, 2] = 3.99990007281303 10^-2; 
    k[4, 2] = 6.49000033736229 10^-2; 
    k[4, 3] = 9.51910100411624 10^-4;

R = 1.987;

Critical Pressure and Temperature and acentric factor for Hydrogen, Methane, Benzene and Toluene

Pc[1] = 190.8; Tc[1] = 59.7; ω[1] = 0.0;

Pc[2] = 673.1; Tc[2] = 343.9; ω[2] = 0.0;

Pc[3] = 714.2; Tc[3] = 1012.7; ω[3] = 0.2116;

Pc[4] = 587.8; Tc[4] = 1069.1; ω[4] = 0.2415;

Feed Pressure and temperature are 485 psia and 100°F. Feed is composed of 0.41% Hydrogen, 5.71% Methane, 70.97% Benzene and 22.91% Toluene. Feed pressure is decreased adiabatically from 485 psia to 165 psia.

P = 165;

z[1] := 0.0041; z[2] = 0.0571; z[3] = 0.7097; z[4] = 0.2291;

z[1] + z[2] + z[3] + z[4]

Adiabatic Flash Calculations using Enthalpy balance with Q=0.

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));

AC[1, 1] = 29.088; AC[2, 1] = -0.192 10^-2; 
    AC[3, 1] = 0.4 10^-5; 
    AC[4, 1] = -0.87 10^-9;
AC[1, 2] = 19.875 ; AC[2, 2] = 5.021 10^-2; 
    AC[3, 2] =  1.268 10^-5; 
    AC[4, 2] =  -11.004 10^-9;
AC[1, 3] = -36.193; AC[2, 3] = 48.444 10^-2; 
    AC[3, 3] = -31.548 10^-5; 
    AC[4, 3] = 77.573 10^-9;
AC[1, 4] = -34.364; AC[2, 4] = 55.887 10^-2; 
    AC[3, 4] = -34.435 10^-5; 
    AC[4, 4] = 80.335 10^-9;

HL1 = (∑_(i = 1)^4) AC[1, 
                i] ((T - 459.67 - 32)/1.8 + 273.15) x[
                i] + (∑_(i = 1)^4) AC[2, 
                i] ((T - 459.67 - 32)/1.8 + 273.15)^2/2 x[
                i] + (∑_(i = 1)^4) AC[3, 
                i] ((T - 459.67 - 32)/1.8 + 273.15)^3/3 x[
                i] + (∑_(i = 1)^4) AC[4, 
                i] ((T - 459.67 - 32)/1.8 + 273.15)^4/4 
               x[i] - ((∑_(i = 1)^4) AC[1, i] (298.15) x[
                    i] + (∑_(i = 1)^4) AC[2, i] (298.15)^2/2 x[
                    i] + (∑_(i = 1)^4) AC[3, i] (298.15)^3/3 x[
                    i] + (∑_(i = 1)^4) AC[4, i] (298.15)^4/4 
                   x[i]);

HL2 = 
        1.987 T ((Subscript[Z, l]) - 1) + 
          1/(2 (Sqrt[2]) 
                  (Subscript[B, l]) 1.987 T/
                    P) Log[((Subscript[Z, l]) + (1 + (Sqrt[2])) 
                      (Subscript[B, l]))/((Subscript[Z, l]) + (1 - (Sqrt[2])) (Subscript[B, l]))] (T D[
                    (Subscript[A, l]) (1.987 T)^2/P, {T, 1}] - 
                (Subscript[A, l]) (1.987 T)^2/P);

HL2 = % /. P -> 165;

HL2 = %/9.486 10^4/453.593;

HL = HL1 + HL2;

HV1 = (∑_(i = 1)^4) AC[1, 
                i] ((T - 459.67 - 32)/1.8 + 273.15) y[
                i] + (∑_(i = 1)^4) AC[2, 
                i] ((T - 459.67 - 32)/1.8 + 273.15)^2/2 y[
                i] + (∑_(i = 1)^4) AC[3, 
                i] ((T - 459.67 - 32)/1.8 + 273.15)^3/3 y[
                i] + (∑_(i = 1)^4) AC[4, 
                i] ((T - 459.67 - 32)/1.8 + 273.15)^4/4 
               y[i] - ((∑_(i = 1)^4) AC[1, i] (298.15) y[
                    i] + (∑_(i = 1)^4) AC[2, i] (298.15)^2/2 y[
                    i] + (∑_(i = 1)^4) AC[3, i] (298.15)^3/3 y[
                    i] + (∑_(i = 1)^4) AC[4, i] (298.15)^4/4 
                   y[i]);

HV2 = 
        1.987 T ((Subscript[Z, v]) - 1) + 
          1/(2 (Sqrt[2]) 
                  (Subscript[B, v]) 1.987 T/
                    P) Log[((Subscript[Z, v]) + (1 + (Sqrt[2])) 
                      (Subscript[B, v]))/((Subscript[Z, v]) + (1 - (Sqrt[2])) (Subscript[B, v]))] (T D[
                    (Subscript[A, v]) (1.987 T)^2/P, {T, 1}] - 
                (Subscript[A, v]) (1.987 T)^2/P);

HV2 = % /. P -> 165;

HV2 = %/9.486 10^4/453.593;

HV = HV1 + HV2;

(100 - 32)/1.8 + 273.15

apf[i_] := 
      0.45724 af[i]  R^2 Tc[i]^2/Pc[i]

Apf[i_] := apf[i] Pf/(R T)^2;

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

Bpf[i_] := Pf bpf[i]/(R T)

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

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

AF = (∑_(i = 1)^4)(∑_(j = 1)^4) z[i] z[j] AFb[i, 
                j];

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

BF = (∑_(i = 1)^4) z[i] Bpf[i];

solZF = 
      Solve[(Subscript[Z, F])^3 + (-1 + BF) (Subscript[Z, F])^2 + 
                  (Subscript[Z, F]) (AF - 3  BF^2 - 2  BF) - AF BF + BF^2 + 
                  BF^3 == 0 /. T ->  100 + 459.67 /. 
            Pf ->  485, (Subscript[Z, F])] // Chop

HF1 = (∑_(i = 1)^4) AC[1, i] (310.928) z[
              i] + (∑_(i = 1)^4) AC[2, i] (310.928)^2/2 z[
              i] + (∑_(i = 1)^4) AC[3, i] (310.928)^3/3 z[
              i] + (∑_(i = 1)^4) AC[4, i] (310.928)^4/4 
             z[i] - ((∑_(i = 1)^4) AC[1, i] (298.15) z[
                  i] + (∑_(i = 1)^4) AC[2, i] (298.15)^2/2 z[
                  i] + (∑_(i = 1)^4) AC[3, i] (298.15)^3/3 z[
                  i] + (∑_(i = 1)^4) AC[4, i] (298.15)^4/4 
                 z[i])

HF2 = 
        1.987 T ((Subscript[Z, F]) - 1) + 
          1/(2 (Sqrt[2]) 
                  BF 1.987 T/
                    Pf) Log[((Subscript[Z, F]) + (1 + (Sqrt[2])) 
                      BF)/((Subscript[Z, F]) + (1 - (Sqrt[2])) BF)] (T D[
                    AF (1.987 T)^2/Pf, {T, 1}] - 
                AF (1.987 T)^2/Pf);

HF2 = % /. Pf -> 485 /. T -> (100 + 459.67) /. 
          solZF[[1]];

HF2 = %/9.486 10^4/453.593

HF = HF1 + HF2

sol = 
      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, 
            HF ==  (Subscript[ϕ, vap]) HV + (1 - (Subscript[ϕ, vap])) HL}, {x[
              1], 0.0004}, {x[2], 0.03}, {x[3], 0.7}, {x[4], 
            0.25}, {(Subscript[ϕ, vap]), 0.01}, {y[1], 0.1}, {y[2], 0.85}, {y[3], 
            0.02}, {y[4], 0.001}, {(Subscript[Z, l]), 0.05}, {(Subscript[Z, v]), 0.95}, {T, 540}, 
          MaxIterations -> 1000] // Chop

Results of the iterative technique of Henley and Seader:

Equilibrium liquid mole fractions:  { 0.000413642, 0.0280863, 0.734359, 0.237141 } ,
Equilibrium vapor mole fractions: { 0.108406, 0.869565, 0.0197101, 0.00231884 } ,
Only 3.44% of the feed is vaporized.

All these results compare very well with the approach using FindRoot presented here.