Residue Curve Map for the mixture 
MTBE/Methanol/isobutene/n-butane at P=1 atm
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Activity Coefficient using Wilson model
Component 1: Isobutene
Component 2: Methanol
Component 3: MTBE
Component 4: n-Butane

G[i_] := Exp[
                1 - Log[(∑_(j = 1)^4) x[j] A[i, 
                          j]] - (∑_(k = 1)^4)(((x[k] A[k, i])/((∑_(j = 1)^4) x[j] A[k, j])))] /. x[3] -> x3 /. 
            x[4] -> x4 /. x[1] -> x1 /. 
        x[2] -> x2

A[i_, j_] := V[j]/V[i] Exp[-a[i, j]/(R T)]

a[1, 1] = 0;
a[1, 2] = 169.9953;
a[1, 3] = -30.2477;
a[1, 4] = 0;
a[2, 1] = 2576.8532;
a[2, 2] = 0;
a[2, 3] = 1483.2478;
a[2, 4] = 2283.8726;
a[3, 1] = 271.5669;
a[3, 2] = -406.3902;
a[3, 3] = 0;
a[3, 4] = 0;
a[4, 1] = 0;
a[4, 2] = 382.3429;
a[4, 3] = 0;
a[4, 4] = 0;

V[1] = 93.33;
V[2] = 44.44;
V[3] = 118.8;
V[4] = 100.39;

Gas Constant and Total Pressure

R = 1.987; P = 1.013 10^5;

Vapor Pressure Using Antoine Equation and Normal Boiling Temperatures

Psat[i_] := Exp[(Aa[i] + Ba[i]/(T + Ca[i]))]

Aa[1] = 20.64559; Aa[2] = 23.49989;
Aa[3] = 20.71616; Aa[4] = 20.56490;
Ba[1] = -2125.74886; Ba[2] = -3643.31362; 
    Ba[3] = -2571.58460; 
    Ba[4] = -2151.97590;
Ca[1] = -33.160; Ca[2] = -33.434; Ca[3] = -48.406; 
    Ca[4] = -34.420;

FindRoot[Psat[1] ==  1.013 10^5, {T, 300}]

FindRoot[Psat[2] ==  1.013 10^5, {T, 300}]

FindRoot[Psat[3] ==  1.013 10^5, {T, 300}]

FindRoot[Psat[4] ==  1.013 10^5, {T, 300}]

Equilibrium Constant for the reaction of isobutene with methanol to produce MTBE with n-butane as an inert component:
	isobutene + C(Subscript[H, 3])OH +n-butane⇌ MTBE + n-butane

Keq = Exp[-10.0982 + 4254.05/T + 0.2667 Log[T]];

Computing the Residue Curve Map takes about 10 minutes on a Pentium 800 MHz PC

t = 0; h = 0.1; X[1, 0] = 0.2; X[2, 0] = 0.2; 
    X[4, 0] = 1 - X[1, 0] - X[2, 0]; x[1, -1] = 0.1; x[2, -1] = 0.2; 
    x[3, -1] = 0.1; x[4, -1] = 0.8;

While[
        t < 750, 
{
sol = 
            Chop[FindRoot[{X[1, 
                      t] == (x1 + x3)/(1 + 
                          x3), 
X[2, 
                      t] ==  (x2 + x3)/(1 + x3), 
                  G[1] G[2] x1 x2 - G[3] x3/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] x1 + 
                      Psat[2] G[2] x2 + Psat[3] G[3] x3 + 
                      Psat[4] G[4] x4 ==  P}, {T, 350}, {x3, 
                  x[3, t - 1]}, {x4, x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, 
                  x[2, t - 1]}, 
                MaxIterations -> 100000]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
x[4, t] = 
            Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, 
              t] = Chop[Psat[1] G[1] x[1, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[2, 
              t] = Chop[Psat[2] G[2] x[2, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[3, 
              t] = Chop[Psat[3] G[3] x[3, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
Y[1, 
              t] = Chop[(y[1, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[2, t] = 
            Chop[(y[2, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[4, t] = 
            1 - Y[1, t] - Y[2, t], 
X[1, t + 1] = 
            Chop[X[1, t] - h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] - h (X[2, t] - Y[2, t])], 
          X[4, t + 1] = 1 - X[1, t + 1] - X[2, t + 1], t++}] // 
      Timing

tble1 = Table[{X[1, t], X[2, t]}, {t, 0, 750}];

Last[tble1]

plt1 = 
      ListPlot[tble1, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[0, 1, 0]]

t = 0; h = 0.1; X[1, 0] = 0.2; X[2, 0] = 0.2; 
    X[4, 0] = 1 - X[1, 0] - X[2, 0]; x[1, -1] = 0.1; x[2, -1] = 0.2; 
    x[3, -1] = 0.1; x[4, -1] = 0.8;

While[
        t < 110, 
{
sol = 
            Chop[FindRoot[{X[1, 
                      t] == (x1 + x3)/(1 + 
                          x3), 
X[2, 
                      t] ==  (x2 + x3)/(1 + x3), 
                  G[1] G[2] x1 x2 - G[3] x3/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] x1 + 
                      Psat[2] G[2] x2 + Psat[3] G[3] x3 + 
                      Psat[4] G[4] x4 ==  P}, {T, 350}, {x3, 
                  x[3, t - 1]}, {x4, x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, 
                  x[2, t - 1]}, 
                MaxIterations -> 100000]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
x[4, t] = 
            Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, 
              t] = Chop[Psat[1] G[1] x[1, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[2, 
              t] = Chop[Psat[2] G[2] x[2, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[3, 
              t] = Chop[Psat[3] G[3] x[3, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
Y[1, 
              t] = Chop[(y[1, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[2, t] = 
            Chop[(y[2, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[4, t] = 
            1 - Y[1, t] - Y[2, t], 
X[1, t + 1] = 
            Chop[X[1, t] + h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] + h (X[2, t] - Y[2, t])], 
          X[4, t + 1] = 1 - X[1, t + 1] - X[2, t + 1], t++}] // 
      Timing

tble2 = Table[{X[1, t], X[2, t]}, {t, 0, 110}];

Last[tble2]

plt2 = 
      ListPlot[tble2, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[0, 1, 0]]

t = 0; h = 0.1; X[1, 0] = 0.2; X[2, 0] = 0.5; 
    X[4, 0] = 1 - X[1, 0] - X[2, 0]; x[1, -1] = 0.1; x[2, -1] = 0.2; 
    x[3, -1] = 0.1; x[4, -1] = 0.8;

While[
        t < 250, 
{
sol = 
            Chop[FindRoot[{X[1, 
                      t] == (x1 + x3)/(1 + 
                          x3), 
X[2, 
                      t] ==  (x2 + x3)/(1 + x3), 
                  G[1] G[2] x1 x2 - G[3] x3/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] x1 + 
                      Psat[2] G[2] x2 + Psat[3] G[3] x3 + 
                      Psat[4] G[4] x4 ==  P}, {T, 350}, {x3, 
                  x[3, t - 1]}, {x4, x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, 
                  x[2, t - 1]}, 
                MaxIterations -> 100000]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
x[4, t] = 
            Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, 
              t] = Chop[Psat[1] G[1] x[1, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[2, 
              t] = Chop[Psat[2] G[2] x[2, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[3, 
              t] = Chop[Psat[3] G[3] x[3, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
Y[1, 
              t] = Chop[(y[1, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[2, t] = 
            Chop[(y[2, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[4, t] = 
            1 - Y[1, t] - Y[2, t], 
X[1, t + 1] = 
            Chop[X[1, t] - h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] - h (X[2, t] - Y[2, t])], 
          X[4, t + 1] = 1 - X[1, t + 1] - X[2, t + 1], t++}] // 
      Timing

tble3 = Table[{X[1, t], X[2, t]}, {t, 0, 250}];

Last[tble3]

plt3 = 
      ListPlot[tble3, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[1, 0, 0]]

t = 0; h = 0.05; X[1, 0] = 0.2; X[2, 0] = 0.5; 
    X[4, 0] = 1 - X[1, 0] - X[2, 0]; x[1, -1] = 0.1; x[2, -1] = 0.2; 
    x[3, -1] = 0.1; x[4, -1] = 0.8;

While[
        t < 60, 
{
sol = 
            Chop[FindRoot[{X[1, 
                      t] == (x1 + x3)/(1 + 
                          x3), 
X[2, 
                      t] ==  (x2 + x3)/(1 + x3), 
                  G[1] G[2] x1 x2 - G[3] x3/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] x1 + 
                      Psat[2] G[2] x2 + Psat[3] G[3] x3 + 
                      Psat[4] G[4] x4 ==  P}, {T, 350}, {x3, 
                  x[3, t - 1]}, {x4, x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, 
                  x[2, t - 1]}, 
                MaxIterations -> 100000]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
x[4, t] = 
            Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, 
              t] = Chop[Psat[1] G[1] x[1, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[2, 
              t] = Chop[Psat[2] G[2] x[2, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[3, 
              t] = Chop[Psat[3] G[3] x[3, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
Y[1, 
              t] = Chop[(y[1, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[2, t] = 
            Chop[(y[2, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[4, t] = 
            1 - Y[1, t] - Y[2, t], 
X[1, t + 1] = 
            Chop[X[1, t] + h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] + h (X[2, t] - Y[2, t])], 
          X[4, t + 1] = 1 - X[1, t + 1] - X[2, t + 1], t++}] // 
      Timing

tble4 = Table[{X[1, t], X[2, t]}, {t, 0, 60}];

Last[tble4]

plt4 = 
      ListPlot[tble4, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[1, 0, 0]]

t = 0; h = 0.1; X[1, 0] = 0.5; X[2, 0] = 0.2; 
    X[4, 0] = 1 - X[1, 0] - X[2, 0]; x[1, -1] = 0.4; x[2, -1] = 0.4; 
    x[3, -1] = 0.1; x[4, -1] = 0.1;

While[
        t < 250, 
{
sol = 
            Chop[FindRoot[{X[1, 
                      t] == (x1 + x3)/(1 + 
                          x3), 
X[2, 
                      t] ==  (x2 + x3)/(1 + x3), 
                  G[1] G[2] x1 x2 - G[3] x3/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] x1 + 
                      Psat[2] G[2] x2 + Psat[3] G[3] x3 + 
                      Psat[4] G[4] x4 ==  P}, {T, 350}, {x3, 
                  x[3, t - 1]}, {x4, x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, 
                  x[2, t - 1]}, 
                MaxIterations -> 100000]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
x[4, t] = 
            Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, 
              t] = Chop[Psat[1] G[1] x[1, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[2, 
              t] = Chop[Psat[2] G[2] x[2, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[3, 
              t] = Chop[Psat[3] G[3] x[3, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
Y[1, 
              t] = Chop[(y[1, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[2, t] = 
            Chop[(y[2, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[4, t] = 
            1 - Y[1, t] - Y[2, t], 
X[1, t + 1] = 
            Chop[X[1, t] - h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] - h (X[2, t] - Y[2, t])], 
          X[4, t + 1] = 1 - X[1, t + 1] - X[2, t + 1], t++}] // 
      Timing

tble5 = Table[{X[1, t], X[2, t]}, {t, 0, 250}];

Last[tble5]

plt5 = 
      ListPlot[tble5, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[0, 0, 1]]

t = 0; h = 0.1; X[1, 0] = 0.5; X[2, 0] = 0.2; 
    X[4, 0] = 1 - X[1, 0] - X[2, 0]; x[1, -1] = 0.4; x[2, -1] = 0.4; 
    x[3, -1] = 0.1; x[4, -1] = 0.1;

While[
        t < 110, 
{
sol = 
            Chop[FindRoot[{X[1, 
                      t] == (x1 + x3)/(1 + 
                          x3), 
X[2, 
                      t] ==  (x2 + x3)/(1 + x3), 
                  G[1] G[2] x1 x2 - G[3] x3/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] x1 + 
                      Psat[2] G[2] x2 + Psat[3] G[3] x3 + 
                      Psat[4] G[4] x4 ==  P}, {T, 350}, {x3, 
                  x[3, t - 1]}, {x4, x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, 
                  x[2, t - 1]}, 
                MaxIterations -> 100000]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
x[4, t] = 
            Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, 
              t] = Chop[Psat[1] G[1] x[1, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[2, 
              t] = Chop[Psat[2] G[2] x[2, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
y[3, 
              t] = Chop[Psat[3] G[3] x[3, t] 1/P] /. 
                      x3 -> x[3, t] /. x4 -> x[4, t] /. x1 -> x[1, t] /. 
                x2 -> x[2, t] /. T -> T[t], 
Y[1, 
              t] = Chop[(y[1, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[2, t] = 
            Chop[(y[2, t] + y[3, t])/(1 + 
                    y[3, t])], 
Y[4, t] = 
            1 - Y[1, t] - Y[2, t], 
X[1, t + 1] = 
            Chop[X[1, t] + h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] + h (X[2, t] - Y[2, t])], 
          X[4, t + 1] = 1 - X[1, t + 1] - X[2, t + 1], t++}] // 
      Timing

tble6 = Table[{X[1, t], X[2, t]}, {t, 0, 110}];

Last[tble6]

plt6 = 
      ListPlot[tble6, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[0, 0, 1]]

Show[plt1, plt2, plt3, plt4, plt5, plt6, 
      PlotRange -> {{0, 1}, {0, 1}}, 
      AspectRatio -> 
        1, 
Epilog -> {Line[{{1, 0}, {0, 
                1}}]}]