Phase Envelope for Ethane-Benzene Mixture 
Authors' Data: Zakia Nasri & 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"]

Off[Unset::"norep"]

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

Phase Envelope for Mixture 1: 
20% Ethane and 80% Benzene

P = 5 14.69;

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

y[1] =. ; y[2] =. ;

T =.

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}, {T, 
            400}, {(Subscript[Z, l]), 0.05}, {(Subscript[Z, v]), 0.9}, MaxIterations -> 1000000] // 
        Chop

Te[1] = T /. sol[1]; Pe[1] = 5 14.69;

i =.

For[i = 2, i < 65, P = (4 + i) 14.69; x[1] = 0.2; 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, 800}, {(Subscript[Z, l]), 
              (Subscript[Z, l]) /. sol[i - 1]}, {(Subscript[Z, v]), (Subscript[Z, v]) /. sol[i - 1]}, 
            MaxIterations -> 1000000] // Chop; Te[i] = T /. sol[i]; 
      Pe[i] = P; 
      Print[i, " ", Te[i], " ", Pe[i], " ", 
        (Subscript[Z, l]) /. sol[i], " ", (Subscript[Z, v]) /. sol[i]]; i++]

tbl1 = 
        Append[Table[{((Te[i] - 459.67) - 32)/1.8 + 273.15, 
              Pe[i]/14.69}, {i, 1, 64}], {((513) - 32)/1.8 + 273.15, 
            968/14.69}];

plt1 = 
      ListPlot[tbl1, PlotJoined -> True, 
        PlotStyle -> RGBColor[1, 0, 0], 
        PlotRange -> {{200, 600}, {0, 100}}, 
        DisplayFunction -> Identity]

P = 5 14.69;

y[1] = 0.2; y[2] = 1 - x[1];

T =.

x[1] =. ; x[2] =. ;

sol[1] = 
      FindRoot[{y[1] ==  K[1] x[1], y[2] ==  K[2] x[2], 
            x[1] + x[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}, {x[1], 0.005}, {x[2], 0.995}, {T, 
            400}, {(Subscript[Z, l]), 0.05}, {(Subscript[Z, v]), 0.9}, MaxIterations -> 1000000] // 
        Chop

Tro[1] = T /. sol[1]; Pr[1] = 5 14.69;

i =.

For[i = 2, i < 62, P = (4 + i) 14.69; y[1] = 0.2; y[2] = 1 - y[1]; 
      sol[i] = FindRoot[{y[1] ==  K[1] x[1], 
              y[2] ==  K[2] x[2], x[1] + x[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}, {x[1], x[1] /. sol[i - 1]}, {x[2], 
              x[2] /. sol[i - 1]}, {T, 800}, {(Subscript[Z, l]), 
              (Subscript[Z, l]) /. sol[i - 1]}, {(Subscript[Z, v]), (Subscript[Z, v]) /. sol[i - 1]}, 
            MaxIterations -> 1000000] // Chop; Tro[i] = T /. sol[i]; 
      Pr[i] = P; 
      Print[i, " ", Tro[i], " ", Pr[i], " ", 
        (Subscript[Z, l]) /. sol[i], " ", (Subscript[Z, v]) /. sol[i]]; i++]

tbl2 = 
        Append[Table[{((Tro[i] - 459.67) - 32)/1.8 + 273.15, 
              Pr[i]/14.69}, {i, 1, 61}], {((513) - 32)/1.8 + 273.15, 
            968/14.69}];

plt2 = 
      ListPlot[tbl2, PlotJoined -> True, 
        PlotStyle -> RGBColor[0, 0, 1], 
        PlotRange -> {{200, 600}, {0, 100}}, 
        DisplayFunction -> Identity]

Show[plt1, plt2, 
      Graphics[{RGBColor[0.75, 0.75, 0], PointSize[0.02], 
          Point[{((513) - 32)/1.8 + 273.15, 968/14.69}]}], 
      DisplayFunction -> $DisplayFunction]

Phase Envelope for Mixture 2: 
70% Ethane and 30% Benzene

P = 5 14.69;

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

y[1] =. ; y[2] =. ;

T =.

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}, {T, 
            400}, {(Subscript[Z, l]), 0.05}, {(Subscript[Z, v]), 0.9}, MaxIterations -> 1000000] // 
        Chop

Te[1] = T /. sol[1]; Pe[1] = 5 14.69;

i =.

For[i = 2, i < 101, P = (4 + i) 14.69; x[1] = 0.7; 
      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; Te[i] = T /. sol[i]; 
      Pe[i] = P; 
      Print[i, " ", Te[i], " ", Pe[i], " ", 
        (Subscript[Z, l]) /. sol[i], " ", (Subscript[Z, v]) /. sol[i]]; i++]

tbl3 = 
        Append[Table[{((Te[i] - 459.67) - 32)/1.8 + 273.15, 
              Pe[i]/14.69}, {i, 1, 
              100}], {((770 - 459.67) - 32)/1.8 + 273.15, 
            1560/14.69}];

plt3 = 
      ListPlot[tbl3, PlotJoined -> True, 
        PlotStyle -> RGBColor[1, 0, 0], 
        PlotRange -> {{200, 600}, {0, 120}}, 
        DisplayFunction -> Identity]

P = 5 14.69;

y[1] = 0.7; y[2] = 1 - x[1];

T =.

x[1] =. ; x[2] =. ;

sol[1] = 
      FindRoot[{y[1] ==  K[1] x[1], y[2] ==  K[2] x[2], 
            x[1] + x[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}, {x[1], 0.005}, {x[2], 0.995}, {T, 
            400}, {(Subscript[Z, l]), 0.05}, {(Subscript[Z, v]), 0.9}, MaxIterations -> 1000000] // 
        Chop

Tro[1] = T /. sol[1]; Pr[1] = 5 14.69;

i =.

For[i = 2, i < 101, P = (4 + i) 14.69; y[1] = 0.7; 
      y[2] = 1 - y[1]; 
      sol[i] = FindRoot[{y[1] ==  K[1] x[1], 
              y[2] ==  K[2] x[2], x[1] + x[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}, {x[1], x[1] /. sol[i - 1]}, {x[2], 
              x[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; Tro[i] = T /. sol[i]; 
      Pr[i] = P; 
      Print[i, " ", Tro[i], " ", Pr[i], " ", 
        (Subscript[Z, l]) /. sol[i], " ", (Subscript[Z, v]) /. sol[i]]; i++]

tbl4 = 
        Append[Table[{((Tro[i] - 459.67) - 32)/1.8 + 273.15, 
              Pr[i]/14.69}, {i, 1, 
              100}], {((770 - 459.67) - 32)/1.8 + 273.15, 
            1560/14.69}];

plt4 = 
      ListPlot[tbl4, PlotJoined -> True, 
        PlotStyle -> RGBColor[0, 0, 1], 
        PlotRange -> {{200, 600}, {0, 120}}, 
        DisplayFunction -> Identity]

Show[plt3, plt4, 
      Graphics[{RGBColor[0.75, 0.75, 0], PointSize[0.02], 
          Point[{((770 - 459.67) - 32)/1.8 + 273.15, 1560/14.69}]}], 
      DisplayFunction -> $DisplayFunction]

Plotting both envelopes together. This mixture present a peculiar phenomenon called retrograde condensation.

Show[plt3, plt4, plt1, plt2, 
      Graphics[{RGBColor[0.75, 0.75, 0], PointSize[0.02], 
          Point[{((770 - 459.67) - 32)/1.8 + 273.15, 1560/14.69}], 
          Point[{((513) - 32)/1.8 + 273.15, 968/14.69}]}], 
      DisplayFunction -> $DisplayFunction, Frame -> True, 
      GridLines -> Automatic]