Dynamics of a binary distillation column
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgment: This computation was inspired from a Matlab program available at the Nonlinear Model Database by John Hedengren (University of Texas). 
For more models, please visit the following link :

http://www.che.utexas.edu/~john/research/models.htm

Off[General::"spell"]

We consider a binary mixture with constant relative 
volatility α=1.6

α = 1.6;

Table[
        y[i_, t_] := 
          x[i, t]*α/(1 + (α - 1)*x[i, t]), {i, 1, 
          32}];

Governing equations in the rectifying section

eq[1] = 
        D[x[1, t], {t, 1}] ==  
          1/hdcond V (y[2, t] - x[1, t]);

eq[i_] := 
      D[x[i, t], {t, 1}] ==  
        1/hdup (L (x[i - 1, t] - x[i, t]) - 
              V (y[i, t] - y[i + 1, t]))

Governing equations in the feed stage

eq[17] = 
        D[x[17, t], {t, 1}] ==  
          hdup (F xF + L x[16, t] - Lp x[17, t] - 
                V (y[17, t] - y[18, t]));

Governing equations in the stripping section

Eq[i_] := 
      D[x[i, t], {t, 1}] ==  
        1/hdup (Lp (x[i - 1, t] - x[i, t]) - 
              V (y[i, t] - y[i + 1, t]))

Eq[32] = 
        D[x[32, t], {t, 1}] ==  
          1/hdreb (Lp x[31, t] - (F - Dist) x[32, t] - 
                V y[32, t]);

Feed flow rate and composition

xF = 0.5;

F = 24/60. ;

Molar hold up in the reboiler, condenser and trays

hdreb = 1.0;

hdup = 0.25;

hdcond = 0.5;

Reflux ratio, Distillate flow rate, Liquid and vapor flow rates
in both sections of the column

Dist = 0.5 F;

L = R*Dist;

Lp = F + L;

V = L + Dist;

R = 3;

Steady state for R=3.0

sol = 
      NDSolve[Join[{eq[1]}, Table[eq[i], {i, 2, 16}], {eq[17]}, 
          Table[Eq[i], {i, 18, 31}], {Eq[32]}, 
          Table[x[i, 0] ==  0.5, {i, 1, 32}]], 
        Table[x[i, t], {i, 1, 32}], {t, 0, 20000}]

Table[x[i, t] /. sol /. t ->  20000, {i, 1, 32}]

Perform a step in the reflux ratio (take R=1.5).

R = 1.5;

sol2 = 
      NDSolve[Join[{eq[1]}, Table[eq[i], {i, 2, 16}], {eq[17]}, 
          Table[Eq[i], {i, 18, 31}], {Eq[32]}, 
          Table[x[i, 0] ==  x[i, t] /. sol /. 
              t ->  20000, {i, 1, 32}]], 
        Table[x[i, t], {i, 1, 32}], {t, 0, 20000}]

Plot[Evaluate[Flatten[Table[x[i, t] /. sol2, {i, 1, 16}]]], {t, 0, 60}, 
      PlotStyle ->  Table[Hue[i/33], {i, 1, 16}], 
      PlotRange -> {0.4, 1}]

Plot[
      Evaluate[Flatten[Table[x[i, t] /. sol2, {i, 17, 32}]]], {t, 0, 60}, 
      PlotStyle ->  Table[Hue[i/33], {i, 17, 32}], 
      PlotRange -> {0, 0.6}]