Dynamics of a continuous distillation column for the separation of a mixture of
Benzene, Toluene and p-Xylene 
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Feed plate location (Subscript[N, f]) and total number of plates (Subscript[N, t]).

(Subscript[N, t]) = 10; (Subscript[N, f]) = 8;

Mass Balance equations for benzene and toluene in condenser

eq[1, 1] = (Subscript[M, 0]) D[(Subscript[X, 1, 1])[t], {t, 1}] == 
          V (Subscript[Y, 1, 2])[t] - (L + Di) (Subscript[X, 1, 1])[t];

eq[2, 1] = (Subscript[M, 0]) D[(Subscript[X, 2, 1])[t], {t, 1}] == 
          V (Subscript[Y, 2, 2])[t] - (L + Di) (Subscript[X, 2, 1])[t];

Mass Balance equations for benzene and toluene in reboiler

eq[1, (Subscript[N, t])] = (Subscript[M, (Subscript[N, t])]) D[(Subscript[X, 1, (Subscript[N, t])])[t], {t, 1}] == 
          L1 (Subscript[X, 1, (Subscript[N, t]) - 1])[t] - 
            V1 (Subscript[Y, 1, (Subscript[N, t])])[t] - (L1 - V1) 
              (Subscript[X, 1, (Subscript[N, t])])[t];

eq[2, (Subscript[N, t])] = (Subscript[M, (Subscript[N, t])]) D[(Subscript[X, 2, (Subscript[N, t])])[t], {t, 1}] == 
          L1 (Subscript[X, 2, (Subscript[N, t]) - 1])[t] - 
            V1 (Subscript[Y, 2, (Subscript[N, t])])[t] - (L1 - V1) 
              (Subscript[X, 2, (Subscript[N, t])])[t];

Mass Balance equations for benzene and toluene in stripping section

For[i = 1 + (Subscript[N, f]), i < (Subscript[N, t]), i += 1, 
        eq[1, i] = 
          H D[(Subscript[X, 1, i])[t], {t, 1}] == 
            L1 ((Subscript[X, 1, i - 1])[t] - (Subscript[X, 1, i])[t]) + 
              V1 ((Subscript[Y, 1, i + 1])[t] - (Subscript[Y, 1, i])[t])];

For[i = 1 + (Subscript[N, f]), i < (Subscript[N, t]), i += 1, 
        eq[2, i] = 
          H D[(Subscript[X, 2, i])[t], {t, 1}] == 
            L1 ((Subscript[X, 2, i - 1])[t] - (Subscript[X, 2, i])[t]) + 
              V1 ((Subscript[Y, 2, i + 1])[t] - (Subscript[Y, 2, i])[t])];

Mass Balance equations for benzene and toluene in rectifying section

For[i = 2, i < (Subscript[N, f]), i += 1, 
        eq[1, i] = 
          H D[(Subscript[X, 1, i])[t], {t, 1}] == 
            L ((Subscript[X, 1, i - 1])[t] - (Subscript[X, 1, i])[t]) + 
              V ((Subscript[Y, 1, i + 1])[t] - (Subscript[Y, 1, i])[t])];

For[i = 2, i < (Subscript[N, f]), i += 1, 
        eq[2, i] = 
          H D[(Subscript[X, 2, i])[t], {t, 1}] == 
            L ((Subscript[X, 2, i - 1])[t] - (Subscript[X, 2, i])[t]) + 
              V ((Subscript[Y, 2, i + 1])[t] - (Subscript[Y, 2, i])[t])];

Mass Balance equations for benzene and toluene in the feed plate

eq[1, (Subscript[N, f])] = 
        H D[(Subscript[X, 1, (Subscript[N, f])])[t], {t, 1}] == 
          F (Subscript[X, 1, f]) + L (Subscript[X, 1, (Subscript[N, f]) - 1])[t] - L1 (Subscript[X, 1, (Subscript[N, f])])[t] + 
            V1 (Subscript[Y, 1, (Subscript[N, f]) + 1])[t] - V (Subscript[Y, 1, (Subscript[N, f])])[t];

eq[2, (Subscript[N, f])] = 
        H D[(Subscript[X, 2, (Subscript[N, f])])[t], {t, 1}] == 
          F (Subscript[X, 2, f]) + L (Subscript[X, 2, (Subscript[N, f]) - 1])[t] - L1 (Subscript[X, 2, (Subscript[N, f])])[t] + 
            V1 (Subscript[Y, 2, (Subscript[N, f]) + 1])[t] - V (Subscript[Y, 2, (Subscript[N, f])])[t];

Reflux ratio

R=5;

Molar flow rates of vapor and liquid phases in the stripping and rectifying sections. Assume saturated liquid feed.

L1=L+F;F=40;V1=150;V=V1;L=V-Di;Di=V/(R+1);

Feed composition

(Subscript[X, 1, f]) = 0.65; (Subscript[X, 2, f]) = 0.25; 
    (Subscript[X, 3, f]) = 1 - (Subscript[X, 1, f]) - (Subscript[X, 2, f]);

Liquid holdup in plates, reboiler and condenser are H, (Subscript[M, 0]) and (Subscript[M, (Subscript[N, t])]).

(Subscript[M, 0]) = 75; H = 10; (Subscript[M, (Subscript[N, t])]) = 150;

Relative volatilities

(Subscript[α, 1]) = 2.34; (Subscript[α, 2]) = 1; (Subscript[α, 3]) = 0.44;

Equilibrium Realtionships. Assume ideal behavior.

For[i = 1, i < (Subscript[N, t]) + 1, i += 1, 
      (Subscript[Y, 1, i])[t] = (Subscript[α, 1]) 
          (Subscript[X, 1, i])[
              t]/((∑_(j = 1)^3) (Subscript[α, j]) 
                  (Subscript[X, j, i])[t])]

For[i = 1, i < (Subscript[N, t]) + 1, i += 1, 
      (Subscript[Y, 2, i])[t] = (Subscript[α, 2]) 
          (Subscript[X, 2, i])[
              t]/((∑_(j = 1)^3) (Subscript[α, j]) 
                  (Subscript[X, j, i])[t])]

Summation rule

For[i = 1, i < (Subscript[N, t]) + 1, 
      i += 1, {Print[i], 
        (Subscript[X, 3, i])[t] = 1 - (Subscript[X, 1, i])[t] - (Subscript[X, 2, i])[t]}]

Equations and initial conditions

equations = 
        Join[Table[eq[1, i], {i, 1, (Subscript[N, t])}], 
          Table[eq[2, i], {i, 1, (Subscript[N, t])}]];

init = 
        Join[Table[(Subscript[X, 1, i])[0] == (Subscript[X, 1, f]), {i, 1, (Subscript[N, t])}], 
          Table[(Subscript[X, 2, i])[0] == (Subscript[X, 2, f]), {i, 1, 
              (Subscript[N, t])}]];

sys=Join[equations,init];

funct = 
        Join[Table[(Subscript[X, 1, i])[t], {i, 1, (Subscript[N, t])}], 
          Table[(Subscript[X, 2, i])[t], {i, 1, (Subscript[N, t])}]];

Solving the set of coupled ODEs

sol=NDSolve[sys,funct,{t,0,200}]

Plotting the steady state values of compositions versus plate number for benzene, toluene and p-xylene

tbl1 = 
      Table[First[(Subscript[X, 1, i])[t] /. sol /. t -> 200], {i, 1, 
          (Subscript[N, t])}]

tbl2 = 
      Table[First[(Subscript[X, 2, i])[t] /. sol /. t -> 200], {i, 1, 
          (Subscript[N, t])}]

tbl3 = 
      Table[First[1 - (Subscript[X, 1, i])[t] - (Subscript[X, 2, i])[t] /. sol /. 
            t -> 200], {i, 1, (Subscript[N, t])}]

plt1 = 
      ListPlot[tbl1, PlotJoined -> True, 
        PlotStyle -> RGBColor[1, 0, 0], 
        PlotRange -> {{0, (Subscript[N, t])}, {0, 1}}]

plt2 = 
      ListPlot[tbl2, PlotJoined -> True, 
        PlotStyle -> RGBColor[0, 1, 0], 
        PlotRange -> {{0, (Subscript[N, t])}, {0, 1}}]

plt3 = 
      ListPlot[tbl3, PlotJoined -> True, 
        PlotStyle -> RGBColor[0, 0, 1], 
        PlotRange -> {{0, (Subscript[N, t])}, {0, 1}}]

Show[plt1,plt2,plt3]

Plotting the compositions at the bottom of the column versus time for benzene, toluene and p-xylene

plt4 = 
      Plot[Evaluate[(Subscript[X, 1, (Subscript[N, t])])[t] /. sol], {t, 0, 200}, 
        PlotStyle -> RGBColor[1, 0, 0], PlotRange -> {0, 1}]

plt5 = 
      Plot[Evaluate[(Subscript[X, 2, (Subscript[N, t])])[t] /. sol], {t, 0, 200}, 
        PlotStyle -> RGBColor[0, 1, 0], PlotRange -> {0, 1}]

plt6 = 
      Plot[Evaluate[1 - (Subscript[X, 1, (Subscript[N, t])])[t] - (Subscript[X, 2, (Subscript[N, t])])[t] /. sol], {t, 0, 
          200}, PlotStyle -> RGBColor[0, 0, 1], 
        PlotRange -> {0, 1}]

Show[plt4,plt5,plt6]

Plotting the compositions at the top of the column versus time for benzene, toluene and p-xylene

plt7 = 
      Plot[Evaluate[(Subscript[X, 1, 1])[t] /. sol], {t, 0, 200}, 
        PlotStyle -> RGBColor[1, 0, 0], 
        PlotRange -> All]

plt8 = 
      Plot[Evaluate[(Subscript[X, 2, 1])[t] /. sol], {t, 0, 200}, 
        PlotStyle -> RGBColor[0, 1, 0], 
        PlotRange -> All]

plt9 = 
      Plot[Evaluate[1 - (Subscript[X, 1, 1])[t] - (Subscript[X, 2, 1])[t] /. sol], {t, 0, 200}, 
        PlotStyle -> RGBColor[0, 0, 1], 
        PlotRange -> All]

Show[plt7,plt8,plt9]