Vapor Liquid Equilibrium Computation using NRTL 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"]

P = 760;

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

Gam[i_]:=Exp[(((∑_(j = 1)^2)( t[j, i] G[j, i] x[j]))/((∑_(l = 1)^2) (G[l, i] x[l])))+(∑_(j = 1)^2) (((G[i, j] x[j])/((∑_(l = 1)^2)( G[l, j] x[l])))(t[i, j]-(((∑_(r = 1)^2)( t[r, j] G[r, j] x[r]))/((∑_(l = 1)^2) (G[l, j] x[l])))))]/.x[1] -> x1/.x [2] -> 1 - x1

Gamp[i_]:=Exp[(((∑_(j = 1)^2)( t[j, i] G[j, i] x[j]))/((∑_(l = 1)^2) (G[l, i] x[l])))+(∑_(j = 1)^2) (((G[i, j] x[j])/((∑_(l = 1)^2)( G[l, j] x[l])))(t[i, j]-(((∑_(r = 1)^2)( t[r, j] G[r, j] x[r]))/((∑_(l = 1)^2) (G[l, j] x[l])))))]/.x[1] -> xp1/.x [2] -> 1 - xp1

G[i_, j_] := Exp[-a[i, j] t[i, j]]

t[i_, j_] := g[i, j]/(R (T + 273.15))

g[1, 1] = 0; g[2, 2] = 0; g[2, 1] = 2794.6660; g[1, 2] = 570.13620; 
    a[1, 1] = 0; a[2, 2] = 0; a[1, 2] = 0.47; a[2, 1] = 0.47;

Vapor Pressure Using Antoine Equation

A1 = 7.36366; B1 = 1305.198; C1 = 173.427;

A2 = 8.07131; B2 = 1730.630; C2 = 233.426;

PS1 = 10^(A1 - B1/(C1 + T));
PS2 = 10^(A2 - B2/(C2 + T));

FindRoot[PS1 ==  P, {T, 100, 120}]

FindRoot[PS2 ==  P, {T, 100, 120}]

Gas Constant and Pressure

R = 1.987; P = 760;

Heterogeneous Azeotrpe

FindRoot[{y1 ==  Gam[1] x1 PS1/P, 
        Gam[1] x1 PS1/P ==  
          Gamp[1] xp1 PS1/P, 
Gamp[1] xp1 PS1 + 
            Gamp[2] (1 - xp1) PS2 ==  P,  
        Gam[1] x1 PS1 + Gam[2] (1 - x1) PS2 ==  
          P}, 
{x1, 0.005, 0.01}, {xp1, 0.35, 0.4}, {y1, 
        0.2, 0.25}, {T, 90, 95}]

Computing VLE data

i = 0;

While[
      i < 101, {T = 
          Evaluate[
            93.2468 + i (117.727 - 93.2468)/100], 
Temp[
            i] = T, 
x1[
            i] = FindRoot[{Gam[1] x1 PS1 + 
                    Gam[2] (1 - x1) PS2 ==  P}, {x1, 0.3, 
                0.95}][[1, 2]], 
xp1[
            i] = FindRoot[{Gam[1] x1 PS1 + 
                    Gam[2] (1 - x1) PS2 ==  P}, {x1, 0.005, 
                0.01}][[1, 2]], 
y1[i] = 
          Gam[1] x1 PS1/P /. x1 -> x1[i], 
yp1[i] = 
          Gam[1] x1 PS1/P /. x1 -> xp1[i], i++ }]

T-xy and xy plots for the 1Butanol-water system

tb1 = Table[{x1[i], Temp[i]}, {i, 1, 100}];

tb2 = Table[{xp1[i], Temp[i]}, {i, 1, 27}];

tb3 = Table[{y1[i], Temp[i]}, {i, 1, 100}];

tb4 = Table[{yp1[i], Temp[i]}, {i, 1, 27}];

tbf1 = Sort[Join[tb1, tb2]];

tbf2 = Sort[Join[tb3, tb4]];

plt1 = 
      ListPlot[tbf1, 
        PlotStyle -> 
          RGBColor[1, 0, 
            0], 
PlotRange -> {{-0.2, 1.2}, {90, 
              120}}, PlotJoined -> True]

plt2 = 
      ListPlot[tbf2, 
        PlotStyle -> 
          RGBColor[0, 0, 
            1], 
PlotRange -> {{-0.2, 1.2}, {90, 
              120}}, PlotJoined -> True]

Show[plt1, plt2]

tb5 = Table[{x1[i], y1[i]}, {i, 1, 100}];

tb6 = Table[{xp1[i], yp1[i]}, {i, 1, 27}];

tbf3 = Sort[Join[tb5, tb6]];

plt3 = 
      ListPlot[tbf3, 
        PlotStyle -> 
          RGBColor[0, 1, 
            0], 
PlotRange -> {{0, 1}, {0, 1}}, 
        PlotJoined -> True, 
        Epilog -> {RGBColor[1, 0, 0], 
            Line[{{0, 0}, {1, 1}}]}]