VLLE Computation using Two-Parameter Margules Model
Case of a Ternary System : 
 Water- Ethanol- Ethyl Acetate
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"]

Vapor pressure and normal boiling points

Psat1 = Exp[A1 + B1/(C1 + T) + D1 Log[T] + E1 T^F1]

Psat2 = Exp[A2 + B2/(C2 + T) + D2 Log[T] + E2 T^F2]

Psat3 = Exp[A3 + B3/(C3 + T) + D3 Log[T] + E3 T^F3]

A1 = 65.92780; B1 = -7227.53; C1 = 0; D1 = -7.17695; 
    E1 = 4.03 10^-6; F1 = 2;

A2 = 86.4860; B2 = -7931.10; C2 = 0; D2 = -10.24980; 
    E2 = 6.39 10^-6; F2 = 2;

A3 = 88.3760; B3 = -7147.86; C3 = 0; D3 = -10.99170; 
    E3 = 8.55 10^-6; F3 = 2;

FindRoot[Psat1 ==  1.01325 10^2, {T, 380}]

FindRoot[Psat2 ==  1.01325 10^2, {T, 380}]

FindRoot[Psat3 ==  1.01325 10^2, {T, 380}]

Two-Parameter Margules Model

gE = (∑_(i = 1)^3)(∑_(j = 1)^3) A[i, j] x[
              i] x[j]^2

% /. A[1, 1] -> 0 /. A[2, 2] ->  0 /. 
      A[3, 3] ->  0

Exp[
      D[(n[1] + n[2] + n[3]) % /. 
              x[1] ->  (n[1]/(n[1] + n[2] + n[3])) /. 
            x[2] ->  (n[2]/(n[1] + n[2] + n[3])) /. 
          x[3] ->  (n[3]/(n[1] + n[2] + n[3])), {n[1], 
          1}]]

% /. (n[1] + n[2] + n[3]) -> A

% /. n[1] ->  x[1] A /. n[2] ->  x[2] A /. 
      n[3] ->  x[3] A

Gam1 = % // Simplify

% /. x[1] ->  xp[2] /. x[2] ->  xp[3] /. 
                  x[3] ->  xp[1] /. A[2, 1] ->  Ap[3, 2] /. 
              A[3, 2] ->  Ap[1, 3] /. A[3, 1] ->  Ap[1, 2] /. 
          A[1, 2] ->  Ap[2, 3] /. A[2, 3] ->  Ap[3, 1] /. 
      A[1, 3] ->  Ap[2, 1]

Gam2 = % /. xp[1] ->  x[1] /. 
                      xp[2] ->  x[2] /. xp[3] ->  x[3] /. 
                  Ap[3, 2] ->  A[3, 2] /. 
                Ap[1, 3] ->  A[1, 3] /. Ap[1, 2] ->  A[1, 2] /. 
            Ap[2, 3] ->  A[2, 3] /. Ap[3, 1] ->  A[3, 1] /. 
        Ap[2, 1] ->  A[2, 1]

Gam1 /. x[1] ->  xp[3] /. 
                    x[2] ->  xp[1] /. x[3] ->  xp[2] /. 
                A[2, 1] ->  Ap[1, 3] /. A[3, 2] ->  Ap[2, 1] /. 
            A[3, 1] ->  Ap[2, 3] /. A[1, 2] ->  Ap[3, 1] /. 
        A[2, 3] ->  Ap[1, 2] /. A[1, 3] ->  Ap[3, 2]

Gam3 = % /. xp[1] ->  x[1] /. 
                      xp[2] ->  x[2] /. xp[3] ->  x[3] /. 
                  Ap[1, 3] ->  A[1, 3] /. 
                Ap[2, 1] ->  A[2, 1] /. Ap[2, 3] ->  A[2, 3] /. 
            Ap[3, 1] ->  A[3, 1] /. Ap[1, 2] ->  A[1, 2] /. 
        Ap[3, 2] ->  A[3, 2]

A[1, 1] = 0.0; A[2, 1] = 1.6266; A[3, 1] =  3.7016;

A[1, 2] = 0.7917; A[2, 2] = 0.0; A[3, 2] = 0.8557;

A[1, 3] = 1.8656 ; A[2, 3] = 0.7476; A[3, 3] = 0.0;

Gam1p = Gam1 /. x[1] ->  xp[1] /. x[2] ->  xp[2] /. 
          x[3] ->  xp[3];

Gam2p = Gam2 /. x[1] ->  xp[1] /. x[2] ->  xp[2] /. 
          x[3] ->  xp[3];

Gam3p = Gam3 /. x[1] ->  xp[1] /. x[2] ->  xp[2] /. 
          x[3] ->  xp[3];

Azeotropes Determination

Homogeneous Binary Azeotrope between Ethanol and Ethyl acetate

FindRoot[{Gam1 x[1] Psat1 + Gam2 x[2] Psat2 + 
              Gam3 x[3] Psat3 ==  1.01325 10^2, 
          x[1] ==  y[1], x[2] ==  y[2], 
          y[1] ==  Gam1 x[1] Psat1/(1.01325 10^2), 
          y[2] ==  Gam2 x[2] Psat2/(1.01325 10^2), 
          x[3] ==  1 - x[1] - x[2]}, {x[1], 0}, {x[2], 0.44}, {y[1], 
          0}, {y[2], 0.4}, {T, 320}, {x[3], 0.7}] // Chop

345.354 - 273.15

Heterogeneous Binary Azeotrope between Water and Ethyl acetate

sol = 
      Minimize[{(Gam1 x[1] Psat1 + Gam2 x[2] Psat2 + 
                    Gam3 x[3] Psat3 -  
                    1.01325 10^2)^2 + 
(Gam1p xp[
                        1] Psat1 + Gam2p xp[2] Psat2 + 
                    Gam3p xp[3] Psat3 -  
                    1.01325 10^2)^2 + 
(x[3] - ( 
                      1 - x[1] - x[2]))^2 + (xp[3] - ( 
                      1 - xp[1] - xp[2]))^2 + (Gam1p xp[1] -  
                    Gam1 x[1])^2 + (Gam2p xp[2] -  
                    Gam2 x[
                        2])^2 + ((xo2 - x[2]) (xp[1] - 
                          x[1]) - (xp[2] - x[2]) (xo1 - 
                          x[1]))^2 + (xo1 -  
                    Gam1 x[1] Psat1/(1.01325 10^2))^2 + (xo2 -  
                    Gam2 x[2] Psat2/(1.01325 10^2))^2 + (xo1 + 
                    xo2 + xo3 - 1)^2, 1 > x[1] >  0, x[2] == 0, 
            1 > xp[1] > 0.9, 320 < T < 370, 0.2 < xo1 < 0.33, 
            xo2 == 0 , 0.6 < xo3 < 0.75, 0 < xp[3] < 1, 0 < x[3] < 1, 
            1 > xp[3] > 0, xp[2] ==  0, 0 < x[2] < 1}, {x[1], x[2], 
            x[3], xp[1], xp[2], xp[3], xo1, T, xo2, xo3}] // Chop

344.162 - 273.15

Cannot detect the Homogeneous Binary Azeotrope between Ethanol and Water and the Ternary azeotrope between Ethanol-Ethyl Acetate-Water

VLLE Determination

xo2 = 0.05; xo3 = 0.55; xo1 = 0.4;

sol1 = 
      Minimize[{(Gam1 x[1] Psat1 + Gam2 x[2] Psat2 + 
                    Gam3 x[3] Psat3 -  
                    1.01325 10^2)^2 + 
(Gam1p xp[
                        1] Psat1 + Gam2p xp[2] Psat2 + 
                    Gam3p xp[3] Psat3 -  
                    1.01325 10^2)^2 + 
(x[3] - ( 
                      1 - x[1] - x[2]))^2 + (xp[3] - ( 
                      1 - xp[1] - xp[2]))^2 + (Gam1p xp[1] -  
                    Gam1 x[1])^2 + (Gam2p xp[2] -  
                    Gam2 x[
                        2])^2 + ((xo2 - x[2]) (xp[1] - 
                          x[1]) - (xp[2] - x[2]) (xo1 - 
                          x[1]))^2 + (xo1 -  
                    Gam1 x[1] Psat1/(1.01325 10^2))^2 + (xo2 -  
                    Gam2 x[2] Psat2/(1.01325 10^2))^2, 
            1 > x[1] >  0, 1 > xp[1] > 0.9, 320 < T < 370, 0 < xp[3] < 1, 
            0 < x[3] < 1, 1 > xp[3] > 0, 0 < x[2] < 1}, {x[1], x[2], x[3], 
            xp[1], xp[2], xp[3], T}] // Chop

xo2 = 0.02; xo3 = 0.58; xo1 = 0.4;

sol2 = 
      Minimize[{(Gam1 x[1] Psat1 + Gam2 x[2] Psat2 + 
                    Gam3 x[3] Psat3 -  
                    1.01325 10^2)^2 + 
(Gam1p xp[
                        1] Psat1 + Gam2p xp[2] Psat2 + 
                    Gam3p xp[3] Psat3 -  
                    1.01325 10^2)^2 + 
(x[3] - ( 
                      1 - x[1] - x[2]))^2 + (xp[3] - ( 
                      1 - xp[1] - xp[2]))^2 + (Gam1p xp[1] -  
                    Gam1 x[1])^2 + (Gam2p xp[2] -  
                    Gam2 x[
                        2])^2 + ((xo2 - x[2]) (xp[1] - 
                          x[1]) - (xp[2] - x[2]) (xo1 - 
                          x[1]))^2 + (xo1 -  
                    Gam1 x[1] Psat1/(1.01325 10^2))^2 + (xo2 -  
                    Gam2 x[2] Psat2/(1.01325 10^2))^2, 
            1 > x[1] >  0, 1 > xp[1] > 0.9, 320 < T < 370, 0 < xp[3] < 1, 
            0 < x[3] < 1, 1 > xp[3] > 0, 0 < x[2] < 1}, {x[1], x[2], x[3], 
            xp[1], xp[2], xp[3], T}] // Chop

xo2 = 0.08; xo3 = 0.52; xo1 = 0.4;

sol3 = 
      Minimize[{(Gam1 x[1] Psat1 + Gam2 x[2] Psat2 + 
                    Gam3 x[3] Psat3 -  
                    1.01325 10^2)^2 + 
(Gam1p xp[
                        1] Psat1 + Gam2p xp[2] Psat2 + 
                    Gam3p xp[3] Psat3 -  
                    1.01325 10^2)^2 + 
(x[3] - ( 
                      1 - x[1] - x[2]))^2 + (xp[3] - ( 
                      1 - xp[1] - xp[2]))^2 + (Gam1p xp[1] -  
                    Gam1 x[1])^2 + (Gam2p xp[2] -  
                    Gam2 x[
                        2])^2 + ((xo2 - x[2]) (xp[1] - 
                          x[1]) - (xp[2] - x[2]) (xo1 - 
                          x[1]))^2 + (xo1 -  
                    Gam1 x[1] Psat1/(1.01325 10^2))^2 + (xo2 -  
                    Gam2 x[2] Psat2/(1.01325 10^2))^2, 
            1 > x[1] >  0, 1 > xp[1] > 0.9, 320 < T < 370, 0 < xp[3] < 1, 
            0 < x[3] < 1, 1 > xp[3] > 0, 0 < x[2] < 1}, {x[1], x[2], x[3], 
            xp[1], xp[2], xp[3], T}] // Chop

xo2 = 0.12; xo3 = 0.48; xo1 = 0.4;

sol4 = 
      Minimize[{(Gam1 x[1] Psat1 + Gam2 x[2] Psat2 + 
                    Gam3 x[3] Psat3 -  
                    1.01325 10^2)^2 + 
(Gam1p xp[
                        1] Psat1 + Gam2p xp[2] Psat2 + 
                    Gam3p xp[3] Psat3 -  
                    1.01325 10^2)^2 + 
(x[3] - ( 
                      1 - x[1] - x[2]))^2 + (xp[3] - ( 
                      1 - xp[1] - xp[2]))^2 + (Gam1p xp[1] -  
                    Gam1 x[1])^2 + (Gam2p xp[2] -  
                    Gam2 x[
                        2])^2 + ((xo2 - x[2]) (xp[1] - 
                          x[1]) - (xp[2] - x[2]) (xo1 - 
                          x[1]))^2 + (xo1 -  
                    Gam1 x[1] Psat1/(1.01325 10^2))^2 + (xo2 -  
                    Gam2 x[2] Psat2/(1.01325 10^2))^2, 
            1 > x[1] >  0, 1 > xp[1] > 0.9, 320 < T < 370, 0 < xp[3] < 1, 
            0 < x[3] < 1, 1 > xp[3] > 0, 0 < x[2] < 1}, {x[1], x[2], x[3], 
            xp[1], xp[2], xp[3], T}] // Chop

xo2 = 0.15; xo3 = 0.45; xo1 = 0.4;

sol5 = 
      Minimize[{(Gam1 x[1] Psat1 + Gam2 x[2] Psat2 + 
                    Gam3 x[3] Psat3 -  
                    1.01325 10^2)^2 + 
(Gam1p xp[
                        1] Psat1 + Gam2p xp[2] Psat2 + 
                    Gam3p xp[3] Psat3 -  
                    1.01325 10^2)^2 + 
(x[3] - ( 
                      1 - x[1] - x[2]))^2 + (xp[3] - ( 
                      1 - xp[1] - xp[2]))^2 + (Gam1p xp[1] -  
                    Gam1 x[1])^2 + (Gam2p xp[2] -  
                    Gam2 x[
                        2])^2 + ((xo2 - x[2]) (xp[1] - 
                          x[1]) - (xp[2] - x[2]) (xo1 - 
                          x[1]))^2 + (xo1 -  
                    Gam1 x[1] Psat1/(1.01325 10^2))^2 + (xo2 -  
                    Gam2 x[2] Psat2/(1.01325 10^2))^2, 
            1 > x[1] >  0, 1 > xp[1] > 0.9, 320 < T < 370, 0 < xp[3] < 1, 
            0 < x[3] < 1, 1 > xp[3] > 0, 0 < x[2] < 1}, {x[1], x[2], x[3], 
            xp[1], xp[2], xp[3], T}] // Chop

plt1 = 
      ListPlot[{{x[3], x[2]} /. sol1[[2]], {x[3], x[2]} /. 
            sol2[[2]], {x[3], x[2]} /. sol3[[2]], {x[3], x[2]} /. 
            sol4[[2]], {x[3], x[2]} /. sol5[[2]]}, 
        PlotStyle -> {PointSize[0.025], RGBColor[1, 0, 0]}, 
        DisplayFunction ->  Identity]

plt3 = 
      ListPlot[{{0.55, Gam2 x[2] Psat2/(1.01325 10^2)} /. 
            sol1[[2]], {0.58, Gam2 x[2] Psat2/(1.01325 10^2)} /. 
            sol2[[2]], {0.52, Gam2 x[2] Psat2/(1.01325 10^2)} /. 
            sol3[[2]], {0.48, Gam2 x[2] Psat2/(1.01325 10^2)} /. 
            sol4[[2]], {0.45, Gam2 x[2] Psat2/(1.01325 10^2)} /. 
            sol5[[2]]}, 
        PlotStyle -> {PointSize[0.025], RGBColor[0, 0, 1]}, 
        DisplayFunction ->  Identity]

plt2 = 
      ListPlot[{{xp[3], xp[2]} /. sol1[[2]], {xp[3], xp[2]} /. 
            sol2[[2]], {xp[3], xp[2]} /. sol3[[2]], {xp[3], xp[2]} /. 
            sol4[[2]], {xp[3], xp[2]} /. sol5[[2]]}, 
        PlotStyle -> {PointSize[0.025], RGBColor[1, 0, 0]}, 
        DisplayFunction ->  Identity]

Show[plt1, plt2, plt3, PlotRange ->  {{0, 1}, {0, 1}}, 
      Epilog ->  {RGBColor[0, 1, 0], 
          Line[{{x[3], x[2]} /. sol1[[2]], {xp[3], xp[2]} /. 
                sol1[[2]]}], 
Line[{{x[3], x[2]} /. 
                sol2[[2]], {xp[3], xp[2]} /. sol2[[2]]}], 
          Line[{{x[3], x[2]} /. sol4[[2]], {xp[3], xp[2]} /. 
                sol4[[2]]}], 
          Line[{{x[3], x[2]} /. sol3[[2]], {xp[3], xp[2]} /. 
                sol3[[2]]}], 
          Line[{{x[3], x[2]} /. sol5[[2]], {xp[3], xp[2]} /. 
                sol5[[2]]}], RGBColor[0, 0, 0], Line[{{1, 0}, {0, 1}}]}, 
      AspectRatio -> 1, 
      DisplayFunction ->  $DisplayFunction]

Boiling Temperatures

{T /. sol2[[2]], T /. sol1[[2]], T /. sol3[[2]], 
        T /. sol4[[2]], T /. sol5[[2]]} - 273.15