Vapor-Liquid Equilibrium Diagram of Heptane-Hexane at 760 mmHg
Using SRK EOS
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Off[General::"spell1"]

Off[FindRoot::"lstol"]

Reduced Pressure and Temperature

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

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

Equilibrium constants using the ϕ-ϕ method and the Soave-Redlich-Kwong Equation of State

Ap[i_] := 0.42747 a[i] Pre[i]/Tre[i]^2

Bp[i_] := 0.08664 Pre[i]/Tre[i]

(Subscript[ϕ, v])[i_] := 
        Exp[((Subscript[Z, v]) - 1) Bp[i]/(Subscript[B, v]) - Log[(Subscript[Z, v]) - (Subscript[B, v])] - 
            (Subscript[A, v])/(Subscript[B, v]) (2 Ap[i]^0.5/(Subscript[A, v])^0.5 - 
                  Bp[i]/(Subscript[B, v])) Log[((Subscript[Z, v]) + (Subscript[B, v]))/(Subscript[Z, v])]];

(Subscript[ϕ, l])[i_] := 
        Exp[((Subscript[Z, l]) - 1) Bp[i]/(Subscript[B, l]) - Log[(Subscript[Z, l]) - (Subscript[B, l])] - 
            (Subscript[A, l])/(Subscript[B, l]) (2 Ap[i]^0.5/(Subscript[A, l])^0.5 - 
                  Bp[i]/(Subscript[B, l])) Log[((Subscript[Z, l]) + (Subscript[B, l]))/(Subscript[Z, l])]];

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

m[i_] := 0.480 + 1.574 ω[i] - 0.176 ω[i]^2

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

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

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

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

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

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

Critical Pressure and Temperature and acentric factor for heptane and hexane.

Pc[1] = 27.37 760/1.01325; 
    Tc[1] = 267.0 + 273.15; ω[1] = 0.34979;

Pc[2] = 30.32 760/1.01325; 
    Tc[2] = 234.7 + 273.15; ω[2] = 0.30070;

Pressure is 760 mmHg.

P = 760;

VLE Calculations at 760 mmHg.

Sol =.

x[1] = 0; x[2] = 1 - x[1];

Sol[1] = 
        FindRoot[{y[1] ==  K[1] x[1], y[2] ==  K[2] x[2], 
              y[1] + y[2] ==  1, 
              (Subscript[Z, l])^3 - (Subscript[Z, l])^2 + (Subscript[Z, l]) ((Subscript[A, l]) - (Subscript[B, l]) - (Subscript[B, l])^2) - 
                  (Subscript[A, l]) (Subscript[B, l]) == 0, 
              (Subscript[Z, v])^3 - (Subscript[Z, v])^2 + (Subscript[Z, v]) ((Subscript[A, v]) - (Subscript[B, v]) - (Subscript[B, v])^2) - 
                  (Subscript[A, v]) (Subscript[B, v]) == 0}, {y[1], 0.999}, {y[2], 0.001}, {T, 
              350}, {(Subscript[Z, l]), 0.01}, {(Subscript[Z, v]), 0.95}, 
            MaxIterations -> 1000000] // Chop;

ye[1] = y[1] /. Sol[1]; Te[1] = T /. Sol[1]; 
    Print[Te[1], "  ", ye[1], " ", (Subscript[Z, l]) /. Sol[1], " ", 
      (Subscript[Z, v]) /. Sol[1]]

342.192 - 273.15

For[i = 2, i < 101, x[1] = 0.01 (i - 1); x[2] = 1 - x[1]; 
      Sol[i] = FindRoot[{y[1] ==  K[1] x[1], 
              y[2] ==  K[2] x[2], y[1] + y[2] ==  1, 
              (Subscript[Z, l])^3 - (Subscript[Z, l])^2 + (Subscript[Z, l]) ((Subscript[A, l]) - (Subscript[B, l]) - (Subscript[B, l])^2) - 
                  (Subscript[A, l]) (Subscript[B, l]) == 0, 
              (Subscript[Z, v])^3 - (Subscript[Z, v])^2 + (Subscript[Z, v]) ((Subscript[A, v]) - (Subscript[B, v]) - (Subscript[B, v])^2) - 
                  (Subscript[A, v]) (Subscript[B, v]) == 0}, {y[1], y[1] /. Sol[i - 1]}, {y[2], 
              y[2] /. Sol[i - 1]}, {T, T /. Sol[i - 1]}, {(Subscript[Z, l]), 
              (Subscript[Z, l]) /. Sol[i - 1]}, {(Subscript[Z, v]), (Subscript[Z, v]) /. Sol[i - 1]}, 
            MaxIterations -> 1000000] // Chop; ye[i] = y[1] /. Sol[i]; 
      Te[i] = T /. Sol[i]; 
      Print[i, " ", Te[i], "  ", ye[i], " ", 
        (Subscript[Z, l]) /. Sol[i], " ", (Subscript[Z, v]) /. Sol[i]]; i++]

371.212 - 273.15

tbl1 = Table[{1 - ye[i], Te[i] - 273.15}, {i, 1, 100}];

tbl2 = 
        Table[{1 - 0.01 (i - 1), Te[i] - 273.15}, {i, 1, 
            100}];

plt1 = 
      ListPlot[tbl1, PlotJoined -> True, PlotRange -> All, 
        PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction -> Identity]

plt2 = 
      ListPlot[tbl2, PlotJoined -> True, PlotRange -> All, 
        PlotStyle -> RGBColor[0, 0, 1], 
        DisplayFunction -> Identity]

Show[plt1, plt2, DisplayFunction -> $DisplayFunction]

VLE Data from DISTIL at 760 mmHg.

TBL1 = {{0.08, 95.10}, {0.22, 89.77}, {0.36, 85.05}, {0.52, 
          80.30}, {0.72, 75.13}}

TBL2 = {{0.1631, 95.10}, {0.3920, 89.77}, {0.5668, 85.05}, {0.7194, 
          80.30}, {0.8611, 75.13}}

plt3 = 
      ListPlot[TBL1, PlotStyle -> {PointSize[0.025], RGBColor[0, 0, 1]}, 
        DisplayFunction -> Identity]

plt4 = 
      ListPlot[TBL2, PlotStyle -> {PointSize[0.025], RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]

Show[plt3, plt4, plt1, plt2, 
      DisplayFunction -> $DisplayFunction]