Residue Curve Map 
for Ternary Systems Using NDSolve
Case of Benzene-Toluene-p-Xylene System
Authors' Data: Housam BINOUS, AÏda Wakad 
and Saïda Ben Achour 
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Off[General::"spell1"]

Gas Constant and Total Pressure

R = 1.987; P = 760;

Vapor Pressure Using Antoine Equation

A3 = 6.87987; B3 = 1196.76; C3 = 219.161;
A2 = 6.95087; B2 = 1342.31; C2 = 219.187;
A1 = 6.99053; B1 = 1453.430; C1 = 215.310;

PS1[t_] := 
        10^(A1 - B1/(C1 + T[t]));
PS2[t_] := 
        10^(A2 - B2/(C2 + T[t]));
PS3[t_] := 10^(A3 - B3/(C3 + T[t]));

Calculation of Residue Curve Map for different initial conditions

y1[t_] = x1[t] PS1[t]/P;
y2[t_] = x2[t] PS2[t]/P;
y3[t_] = x3[t] PS3[t]/P;

x1[0] = 0.2; x2[0] = 0.2; x3[0] = 0.6;

FindRoot[
      P == 
        Evaluate[x1[0] PS1[0]  + x2[0] PS2[0] + x3[0] PS3[0] ], {T[0], 
        100}]

x1[0] =. ; x2[0] =. ; x3[0] =. ;

sol1 = 
      NDSolve[{D[x1[t], {t, 1}] ==  y1[t] - x1[t], 
          D[x2[t], {t, 1}] ==  y2[t] - x2[t], 
          D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t] + 
                x3[t] PS3[t]  - P) == 0 , 
          T[0] == 91.47506015645139` , x1[0] ==  0.2, 
          x2[0] ==  0.2, x3[0] ==  0.6}, {x1, x2, x3, T}, {t, 
          0, 40}]

plt1 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol1], {t, 0, 40}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.5], 
        DisplayFunction ->  Identity]

sol2 = 
      NDSolve[{-D[x1[t], {t, 1}] ==  
            y1[t] - x1[t], -D[x2[t], {t, 1}] ==  
            y2[t] - x2[t], -D[x3[t], {t, 1}] ==  
            y3[t] - x3[t],  (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t]  - P) ==  0, 
          T[0] == 91.47506015645139` , x1[0] ==  0.2, 
          x2[0] ==  0.2, x3[0] ==  0.6}, {x1, x2, x3, T}, {t, 
          0, 10}]

plt2 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol2], {t, 0, 10}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.5], 
        DisplayFunction ->  Identity]

x1[0] = 0.2; x2[0] = 0.3; x3[0] = 0.5;

FindRoot[
      P == 
        Evaluate[x1[0] PS1[0]  + x2[0] PS2[0] + x3[0] PS3[0]], {T[0], 
        100}]

x1[0] =. ; x2[0] =. ; x3[0] =. ;

sol3 = 
      NDSolve[{-D[x1[t], {t, 1}] ==  
            y1[t] - x1[t], -D[x2[t], {t, 1}] ==  
            y2[t] - x2[t], -D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t]  - P) ==  0, 
          T[0] == 94.47049809606979`, x1[0] ==  0.2, 
          x2[0] ==  0.3, x3[0] ==  0.5}, {x1, x2, x3, T}, {t, 
          0, 20}]

plt3 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol3], {t, 0, 20}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.55], 
        DisplayFunction ->  Identity]

sol4 = 
      NDSolve[{D[x1[t], {t, 1}] ==  y1[t] - x1[t], 
          D[x2[t], {t, 1}] ==  y2[t] - x2[t], 
          D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t]  - P) == 0, 
          T[0] == 94.47049809606979` , x1[0] ==  0.2, 
          x2[0] ==  0.3, x3[0] ==  0.5}, {x1, x2, x3, T}, {t, 
          0, 30}]

plt4 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol4], {t, 0, 30}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.55], 
        DisplayFunction ->  Identity]

x1[0] = 0.2; x2[0] = 0.4; x3[0] = 0.4;

FindRoot[
      P == 
        Evaluate[x1[0] PS1[0] + x2[0] PS2[0]  + x3[0] PS3[0] ], {T[0], 
        100}]

x1[0] =. ; x2[0] =. ; x3[0] =. ;

sol5 = 
      NDSolve[{-D[x1[t], {t, 1}] ==  
            y1[t] - x1[t], -D[x2[t], {t, 1}] ==  
            y2[t] - x2[t], -D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t]  - P) ==  0 , 
          T[0] == 97.74716819532348`, x1[0] ==  0.2, 
          x2[0] ==  0.4, x3[0] ==  0.4}, {x1, x2, x3, T}, {t, 
          0, 20}]

plt5 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol5], {t, 0, 20}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.6], 
        DisplayFunction ->  Identity]

sol6 = 
      NDSolve[{D[x1[t], {t, 1}] ==  y1[t] - x1[t], 
          D[x2[t], {t, 1}] ==  y2[t] - x2[t], 
          D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t] - P) == 0, 
          T[0] == 97.74716819532348`, x1[0] ==  0.2, 
          x2[0] ==  0.4, x3[0] ==  0.4}, {x1, x2, x3, T}, {t, 
          0, 30}]

plt6 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol6], {t, 0, 30}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.6], 
        DisplayFunction ->  Identity]

x1[0] = 0.2; x2[0] = 0.5; x3[0] = 0.3;

FindRoot[
      P == 
        Evaluate[x1[0] PS1[0]  + x2[0] PS2[0] + x3[0] PS3[0] ], {T[0], 
        100}]

x1[0] =. ; x2[0] =. ; x3[0] =. ;

sol7 = 
      NDSolve[{D[x1[t], {t, 1}] ==  y1[t] - x1[t], 
          D[x2[t], {t, 1}] ==  y2[t] - x2[t], 
          D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t] + 
                x3[t] PS3[t]  - P) == 0 , 
          T[0] == 101.35081806971178` , x1[0] ==  0.2, 
          x2[0] ==  0.5, x3[0] ==  0.3}, {x1, x2, x3, T}, {t, 
          0, 40}]

plt7 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol7], {t, 0, 40}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.65], 
        DisplayFunction ->  Identity]

sol8 = 
      NDSolve[{-D[x1[t], {t, 1}] ==  
            y1[t] - x1[t], -D[x2[t], {t, 1}] ==  
            y2[t] - x2[t], -D[x3[t], {t, 1}] ==  
            y3[t] - x3[t],  (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t]  - P) ==  0, 
          T[0] == 101.35081806971178` , x1[0] ==  0.2, 
          x2[0] ==  0.5, x3[0] ==  0.3}, {x1, x2, x3, T}, {t, 
          0, 10}]

plt8 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol8], {t, 0, 10}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.65], 
        DisplayFunction ->  Identity]

x1[0] = 0.2; x2[0] = 0.6; x3[0] = 0.2;

FindRoot[
      P == 
        Evaluate[x1[0] PS1[0]  + x2[0] PS2[0] + x3[0] PS3[0]], {T[0], 
        100}]

x1[0] =. ; x2[0] =. ; x3[0] =. ;

sol9 = 
      NDSolve[{-D[x1[t], {t, 1}] ==  
            y1[t] - x1[t], -D[x2[t], {t, 1}] ==  
            y2[t] - x2[t], -D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t]  - P) ==  0, 
          T[0] == 105.33716156339065`, x1[0] ==  0.2, 
          x2[0] ==  0.6, x3[0] ==  0.2}, {x1, x2, x3, T}, {t, 
          0, 20}]

plt9 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol9], {t, 0, 20}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.7], 
        DisplayFunction ->  Identity]

sol10 = 
      NDSolve[{D[x1[t], {t, 1}] ==  y1[t] - x1[t], 
          D[x2[t], {t, 1}] ==  y2[t] - x2[t], 
          D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t]  - P) == 0, 
          T[0] == 105.33716156339065` , x1[0] ==  0.2, 
          x2[0] ==  0.6, x3[0] ==  0.2}, {x1, x2, x3, T}, {t, 
          0, 30}]

plt10 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol10], {t, 0, 30}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.7], 
        DisplayFunction ->  Identity]

x1[0] = 0.2; x2[0] = 0.7; x3[0] = 0.1;

FindRoot[
      P == 
        Evaluate[x1[0] PS1[0] + x2[0] PS2[0]  + x3[0] PS3[0] ], {T[0], 
        100}]

x1[0] =. ; x2[0] =. ; x3[0] =. ;

sol11 = 
      NDSolve[{-D[x1[t], {t, 1}] ==  
            y1[t] - x1[t], -D[x2[t], {t, 1}] ==  
            y2[t] - x2[t], -D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t]  - P) ==  0 , 
          T[0] == 109.77427562442516`, x1[0] ==  0.2, 
          x2[0] ==  0.7, x3[0] ==  0.1}, {x1, x2, x3, T}, {t, 
          0, 20}]

plt11 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol11], {t, 0, 20}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.75], 
        DisplayFunction ->  Identity]

sol12 = 
      NDSolve[{D[x1[t], {t, 1}] ==  y1[t] - x1[t], 
          D[x2[t], {t, 1}] ==  y2[t] - x2[t], 
          D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t] - P) == 0, 
          T[0] == 109.77427562442516`, x1[0] ==  0.2, 
          x2[0] ==  0.7, x3[0] ==  0.1}, {x1, x2, x3, T}, {t, 
          0, 30}]

plt12 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol12], {t, 0, 30}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.75], 
        DisplayFunction ->  Identity]

x1[0] = 0.1; x2[0] = 0.675; x3[0] = 0.225;

FindRoot[
      P == 
        Evaluate[x1[0] PS1[0] + x2[0] PS2[0]  + x3[0] PS3[0] ], {T[0], 
        100}]

x1[0] =. ; x2[0] =. ; x3[0] =. ;

sol13 = 
      NDSolve[{-D[x1[t], {t, 1}] ==  
            y1[t] - x1[t], -D[x2[t], {t, 1}] ==  
            y2[t] - x2[t], -D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t]  - P) ==  0 , 
          T[0] == 102.69203561815846`, x1[0] ==  0.1, 
          x2[0] ==  0.675, x3[0] ==  0.225}, {x1, x2, x3, 
          T}, {t, 0, 20}]

plt13 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol13], {t, 0, 20}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.725], 
        DisplayFunction ->  Identity]

sol14 = 
      NDSolve[{D[x1[t], {t, 1}] ==  y1[t] - x1[t], 
          D[x2[t], {t, 1}] ==  y2[t] - x2[t], 
          D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t] - P) == 0, 
          T[0] == 102.69203561815846`, x1[0] ==  0.1, 
          x2[0] ==  0.675, x3[0] ==  0.225}, {x1, x2, x3, 
          T}, {t, 0, 30}]

plt14 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol14], {t, 0, 30}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.725], 
        DisplayFunction ->  Identity]

x1[0] = 0.15; x2[0] = 0.7; x3[0] = 0.15;

FindRoot[
      P == 
        Evaluate[x1[0] PS1[0] + x2[0] PS2[0]  + x3[0] PS3[0] ], {T[0], 
        100}]

x1[0] =. ; x2[0] =. ; x3[0] =. ;

sol15 = 
      NDSolve[{-D[x1[t], {t, 1}] ==  
            y1[t] - x1[t], -D[x2[t], {t, 1}] ==  
            y2[t] - x2[t], -D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t]  - P) ==  0 , 
          T[0] == 106.60327639229261`, x1[0] ==  0.15, 
          x2[0] ==  0.7, x3[0] ==  0.15}, {x1, x2, x3, T}, {t, 
          0, 20}]

plt15 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol15], {t, 0, 20}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.735], 
        DisplayFunction ->  Identity]

sol16 = 
      NDSolve[{D[x1[t], {t, 1}] ==  y1[t] - x1[t], 
          D[x2[t], {t, 1}] ==  y2[t] - x2[t], 
          D[x3[t], {t, 1}] ==  
            y3[t] - x3[t], (x1[t] PS1[t]  + x2[t] PS2[t]  + 
                x3[t] PS3[t] - P) == 0, 
          T[0] == 106.60327639229261`, x1[0] ==  0.15, 
          x2[0] ==  0.7, x3[0] ==  0.15}, {x1, x2, x3, T}, {t, 
          0, 30}]

plt16 = 
      ParametricPlot[Evaluate[{x1[t], x2[t]} /. sol16], {t, 0, 30}, 
        AspectRatio -> 1, PlotStyle -> RGBColor[0, 1, 0.725], 
        DisplayFunction ->  Identity]

PLT = 
      Show[plt1, plt2, plt3, plt4, plt5, plt6, plt7, plt8, plt9, plt10, 
        plt11, plt12, plt13, plt14, plt15, plt16, 
        Graphics[Line[{{0, 1}, {1, 0}}]], PlotRange -> All, 
        DisplayFunction -> Identity]

PLT1 = 
      Show[PLT, 
        Graphics[{PointSize[0.025], RGBColor[1, 0, 1], 
            Point[Flatten[{x1[t], x2[t]} /. sol13 /. t ->  1.2]]}], 
        DisplayFunction -> Identity]

Calculation of Pinch Point Curve for a given distillate composition. Pinch Point curve is the red curve in plot.

Flatten[{x1[t], x2[t]} /. sol13 /. t ->  1.2]

lst10 = 
        Table[Flatten[{x1[t], x2[t]} /. sol10], {t, 0, 20, 
            0.1}] ;

f10 = Interpolation[lst10]

FindRoot[
      D[f10[x1], {x1, 
            1}] ==  (First[(x2[t] /. sol13 /. t -> 1.2)] - 
              f10[x1])/(First[(x1[t] /. sol13 /. t -> 1.2)] - 
              x1), {x1, 0.0001}]

solt1 = 
      FindRoot[(x1[t] /. sol10) ==  x1 /. %, {t, 0.1, 
          10}]

lst7 = 
        Table[Flatten[{x1[t], x2[t]} /. sol7], {t, 0, 10, 
            0.1}] ;

f7 = Interpolation[lst7]

FindRoot[
      D[f7[x1], {x1, 
            1}] ==  (First[(x2[t] /. sol13 /. t -> 1.2)] - 
              f7[x1])/(First[(x1[t] /. sol13 /. t -> 1.2)] - 
              x1), {x1, 0.001}]

solt2 = 
      FindRoot[(x1[t] /. sol7) ==  x1 /. %, {t, 0.1, 
          10}]

lst4 = 
        Table[Flatten[{x1[t], x2[t]} /. sol4], {t, 0, 10, 
            0.1}] ;

f4 = Interpolation[lst4]

FindRoot[
      D[f4[x1], {x1, 
            1}] ==  (First[(x2[t] /. sol13 /. t -> 1.2)] - 
              f4[x1])/(First[(x1[t] /. sol13 /. t -> 1.2)] - 
              x1), {x1, 0.001}]

solt3 = 
      FindRoot[(x1[t] /. sol4) ==  x1 /. %, {t, 0.1, 
          10}]

lst1 = 
        Table[Flatten[{x1[t], x2[t]} /. sol1], {t, 0, 40, 
            0.01}] ;

f1 = Interpolation[lst1]

FindRoot[
      D[f1[x1], {x1, 
            1}] ==  (First[(x2[t] /. sol13 /. t -> 1.2)] - 
              f1[x1])/(First[(x1[t] /. sol13 /. t -> 1.2)] - 
              x1), {x1, 0.000005}]

solt4 = 
      FindRoot[(x1[t] /. sol1) == x1 /. %, {t, 0.1, 40}]

lst6 = 
        Table[Flatten[{x1[t], x2[t]} /. sol6], {t, 0, 10, 
            0.001}] ;

f6 = Interpolation[lst6]

FindRoot[
      D[f6[x1], {x1, 
            1}] ==  (First[(x2[t] /. sol13 /. t -> 1.2)] - 
              f6[x1])/(First[(x1[t] /. sol13 /. t -> 1.2)] - 
              x1), {x1, 0.0001}]

solt5 = 
      FindRoot[(x1[t] /. sol6) == x1 /. %, {t, 0.1, 10}]

lst12 = 
        Table[Flatten[{x1[t], x2[t]} /. sol12], {t, 0, 10, 
            0.001}] ;

f12 = Interpolation[lst12]

FindRoot[
      D[f12[x1], {x1, 
            1}] ==  (First[(x2[t] /. sol13 /. t -> 1.2)] - 
              f12[x1])/(First[(x1[t] /. sol13 /. t -> 1.2)] - 
              x1), {x1, 0.0001}]

solt6 = 
      FindRoot[(x1[t] /. sol12) == x1 /. %, {t, 0.1, 
          10}]

lst15 = 
        Table[Flatten[{x1[t], x2[t]} /. sol15], {t, 0, 20, 
            0.001}] ;

f15 = Interpolation[lst15]

FindRoot[
      D[f15[x1], {x1, 
            1}] ==  (First[(x2[t] /. sol13 /. t -> 1.2)] - 
              f15[x1])/(First[(x1[t] /. sol13 /. t -> 1.2)] - 
              x1), {x1, 0.2}]

solt7 = 
      FindRoot[(x1[t] /. sol15) == x1 /. %, {t, 0.1, 
          10}]

PLT2 = 
      Show[PLT1, 
        Graphics[{PointSize[0.025], RGBColor[0, 1, 1], 
            Point[Flatten[{x1[t], x2[t]} /. sol10 /. t ->  t /. 
                  solt1]], 
            Point[Flatten[{x1[t], x2[t]} /. sol7 /. t -> t /.  
                  solt2]], 
            Point[Flatten[{x1[t], x2[t]} /. sol4 /. t ->  t /. 
                  solt3]], 
            Point[Flatten[{x1[t], x2[t]} /. sol1 /. t ->  t /. 
                  solt4]], 
            Point[Flatten[{x1[t], x2[t]} /. sol6 /. t ->  t /. 
                  solt5]], 
            Point[Flatten[{x1[t], x2[t]} /. sol12 /. t ->  t /. 
                  solt6]], 
            Point[Flatten[{x1[t], x2[t]} /. sol15 /. t ->  t /. 
                  solt7]]}], DisplayFunction -> Identity]

tbl = {{0, 0}, 
        Flatten[{x1[t], x2[t]} /. sol10 /. t ->  t /. solt1], 
        Flatten[{x1[t], x2[t]} /. sol7 /. t -> t /.  solt2], 
        Flatten[{x1[t], x2[t]} /. sol1 /. t ->  t /. solt4], 
        Flatten[{x1[t], x2[t]} /. sol4 /. t ->  t /. solt3], 
        Flatten[{x1[t], x2[t]} /. sol6 /. t ->  t /. solt5], 
        Flatten[{x1[t], x2[t]} /. sol12 /. t ->  t /. solt6], 
        Flatten[{x1[t], x2[t]} /. sol15 /. t ->  t /. solt7], 
        Flatten[{x1[t], x2[t]} /. sol13 /. t ->  1.2]}

f = Interpolation[tbl, InterpolationOrder ->  1]

plt17b = 
      Plot[f[x], {x, 0, 0.203409}, PlotStyle ->  RGBColor[1, 0, 0], 
        DisplayFunction ->  Identity]

Show[PLT2, plt17b, 
      DisplayFunction ->  $DisplayFunction]

Calculation of Operating Leaf for a given distillate composition. Operating leaf is the gray region.

<< Graphics`FilledPlot`

lst11 = 
        Join[Table[Flatten[{x1[t], x2[t]} /. sol13], {t, 0, 10, 0.01}], 
          Table[Flatten[{x1[t], x2[t]} /. sol14], {t, 0.01, 10, 
              0.01}]] ;

f11 = Interpolation[lst11]

plt18 = 
      FilledPlot[{f[x], f11[x]}, {x, 0.0001, 0.2}, 
        Fills -> GrayLevel[0.75], 
        DisplayFunction ->  Identity]

Show[PLT1, plt18, PlotRange -> All, AspectRatio ->  1.1, 
      DisplayFunction ->  $DisplayFunction]