Vapor Liquid Equilibrium Computation using UNIQUAC Model: Case of Partially Miscible Binary System 1-Butanol - Water 
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)]

Psat2 = Exp[A2 + B2/(C2 + T)]

A1 = 23.2256; B1 = -3835.18; C1 = -45.343;

A2 = 21.9783; B2 = -3080.66; C2 = -96.150;

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

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

Activity coefficients using UNIQUAC model
Component 1 = Water
Component 2 = 1-Butanol

tau[j_, i_] := Exp[-(U[j, i] - U[i, i])/(R T)]

U[1, 1] = 0.0; U[1, 2] = 581.1471;

U[2, 1] = 68.0083; U[2, 2] = 0.0;

R = 1.987;

Gam1 = Exp[Log[Gam1C] + Log[Gam1R]];

Gam2 = Exp[Log[Gam2C] + Log[Gam2R]];

Gam1C = 
        Exp[Log[phi[1]/x[1]] + z/2 q[1] Log[v[1]/phi[1]] + l[1] - 
            phi[1]/x[1] (∑_(j = 1)^2) x[j] l[j]];

Gam2C = 
        Exp[Log[phi[2]/x[2]] + z/2 q[2] Log[v[2]/phi[2]] + l[2] - 
            phi[2]/x[2] (∑_(j = 1)^2) x[j] l[j]];

r[1] = 0.92; r[2] = 3.4543;

q[1] = 1.4; q[2] = 3.0552;

l[i_] := z/2 (r[i] - q[i]) - (r[i] - 1)

z = 10;

Gam1R = 
        Exp[q[1] (1 - 
                Log[(∑_(j = 1)^2) v[j] tau[j, 
                        1]] - (∑_(j = 1)^2)(( v[j] tau[1, j])/((∑_(k = 1)^2) v[k] tau[k, j])))];

Gam2R = 
        Exp[q[2] (1 - 
                Log[(∑_(j = 1)^2) v[j] tau[j, 
                        2]] - (∑_(j = 1)^2)(( v[j] tau[2, j])/((∑_(k = 1)^2) v[k] tau[k, j])))];

phi[i_] := r[i] x[i]/(∑_(j = 1)^2) r[j] x[j]

v[i_] := q[i] x[i]/(∑_(j = 1)^2) q[j] x[j]

Computing VLE data

For[i = 1, 
      i < 100, {x[1] = i 0.01, 
        x[2] = 1 - x[1], 
sol = 
          FindRoot[
            1.01325 10^5 ==  
              Psat1 x[1] Gam1 + Psat2 x[2] Gam2, {T, 
              380}], 
xe1[i] = x[1], xe2[i] = 1 - x[1], 
        Te[i] = sol[[1, 2]], 
        ye1[i] = Psat1 x[1] Gam1/(1.01325 10^5) /. sol, Print[i], 
        i += 1}]

Plotting VLE data. Shape of equilibrium diagram strongly suggests that system exhibits liquid-liquid behavior.

plt1 = 
      ListPlot[Table[{1 - xe1[i], 1 - ye1[i]}, {i, 1, 99}], 
        PlotJoined ->  True, PlotStyle -> RGBColor[1, 0, 0], 
        PlotRange -> {{0, 1}, {0, 1}}]

Determination of liquid-liquid behavior

x[1] =.

x[2] =.

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

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

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

xp[2] = 1 - xp[1];

min = 
      Minimize[{(x[1] Gam1  -  
                  xp[1] Gam1p)^2 + (1.01325 10^5 -  (Psat1 x[
                          1] Gam1 + 
                      Psat2 x[
                          2] Gam2))^2 + (1.01325 10^5 -  (Psat1 xp[1] Gam1p + Psat2 xp[2] Gam2p))^2, x[1] < xp[1], 
          0.9 < xp[1] < 0.999, 0.4 < x[1] < 0.6, T > 300}, {x[1], xp[1], 
          T}]

Psat1 x[1] Gam1/(1.01325 10^5) /. min[[2]]

Show[plt1, 
      Graphics[{RGBColor[0, 0, 1], 
          Line[{{1 - 0.46507, 1 - 0.759129}, {1 - 0.98160, 
                1 - 0.759129}}]}]]

Plotting the VLE diagram obtained from UNIQUAC model

For[i = 1, i < 100, If[xe1[i] > 0.46507, {Print[i], Break[]}], 
      i += 1]

For[i = 1, i < 100, If[xe1[i] > 0.9816, {Print[i], Break[]}], 
      i += 1]

plt2b = 
      ListPlot[Join[{{1, 1}}, Table[{1 - xe1[i], 1 - ye1[i]}, {i, 1, 47}]], 
        PlotJoined ->  True, PlotStyle -> RGBColor[1, 0, 0], 
        PlotRange -> {{0, 1}, {0, 1}}, 
        DisplayFunction ->  Identity]

plt2c = 
      ListPlot[Join[Table[{1 - xe1[i], 1 - ye1[i]}, {i, 98, 99}], {{0, 0}}], 
        PlotJoined ->  True, PlotStyle -> RGBColor[1, 0, 0], 
        PlotRange -> {{0, 1}, {0, 1}}, 
        DisplayFunction ->  Identity]

Show[plt2c, plt2b, 
      Graphics[{RGBColor[0, 0, 1], 
          Line[{{1 - 0.46507, 1 - 0.759129}, {1 - 0.98160, 1 - 0.759129}}], 
          RGBColor[0, 1, 0], Line[{{0, 0}, {1, 1}}]}], Frame ->  True, 
      DisplayFunction ->  $DisplayFunction]