Vapor-Liquid Equilibrium Diagram of Ethane-Benzene at θ=25°C and 175°C
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 Thermodynamique, Application au Génie Chimique et à l'Industrie Pétrolière by Jean Vidal, Technip, 1997.

Off[General::"spell1"]

Off[FindRoot::"lstol"]

Reduced Pressure and Temperature

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

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

Equilibrium constants using the ϕ-ϕ methode 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 Ethane and Benzene

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

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

Temperature is 77 °F.

25 1.8 + 32

T = 25 1.8 + 32 + 459.67;

VLE Calculations at 77°F.

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}, {P, 
              150}, {(Subscript[Z, l]), 0.1}, {(Subscript[Z, v]), 0.9}, MaxIterations -> 1000000] // 
          Chop;

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

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]}, {P, P /. 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]; 
      Pe[i] = P /. Sol[i]; 
      Print[i, " ", Pe[i], "  ", ye[i], " ", 
        (Subscript[Z, l]) /. Sol[i], " ", (Subscript[Z, v]) /. Sol[i]]; i++]

tbl1 = Table[{ye[i], Pe[i]/14.69}, {i, 1, 100}];

tbl2 = 
        Table[{0.01 (i - 1), Pe[i]/14.69}, {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]

Temperature is 347 °F.

175 1.8 + 32

T = 175 1.8 + 32 + 459.67;

VLE Calculations at 347°F.

Sol =.

x[1] = 0.0001; 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.005}, {y[2], 0.995}, {P, 
            155}, {(Subscript[Z, l]), 0.05}, {(Subscript[Z, v]), 0.9}, MaxIterations -> 1000000] // 
        Chop

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

For[i = 2, i < 68, 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]}, {P, P /. 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]; 
      Pe[i] = P /. Sol[i]; 
      Print[i, " ", Pe[i], "  ", ye[i], " ", 
        (Subscript[Z, l]) /. Sol[i], " ", (Subscript[Z, v]) /. Sol[i]]; i++]

tbl1 = Table[{ye[i], Pe[i]/14.69}, {i, 1, 67}];

tbl2 = 
        Table[{0.01 (i - 1), Pe[i]/14.69}, {i, 1, 67}];

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]