Dynamics of a binary distillation column:
Case of a Cyclohexane-Heptane mixture
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"]

Vapor Pressures and activity coefficients using Wilson's model

Table[
        PS1[i_, t_] := 
          Exp[A1 + B1/T[i, t] + C1 Log[T[i, t]] + D1 T[i, t]^E1], {i, 1, 
          32}];

Table[
        PS2[i_, t_] := 
          Exp[A2 + B2/T[i, t] + C2 Log[T[i, t]] + D2 T[i, t]^E2], {i, 1, 
          32}];

A1 = 5.1087 10^1; A2 = 8.7829 10^1;
B1 = -5.2264 10^3; 
    B2 = -6.9964 10^3;
C1 = -4.2278; C2 = -9.8802;
D1 = 9.7554 10^-18; 
    D2 = 7.2099 10^-6;
E1 = 6.0;E2 = 
      2.0;

Table[
        G1[i_, t_] := 
          Exp[-Log[
                  x[i, t] + A12 (1 - x[i, t])] + (1 - 
                    x[i, t]) (A12/(x[i, t] + 
                          A12 (1 - x[i, t])) - 
                    A21/(A21 x[i, t] + 1 - x[i, t]))], {i, 1, 
          32}];

Table[
        G2[i_, t_] := 
          Exp[-Log[
                  1 - x[i, t] + A21 (x[i, t])] + (x[i, 
                    t]) (A21/(1 - x[i, t] + A21 (x[i, t])) - 
                    A12/(A12 (1 - x[i, t]) + x[i, t]))], {i, 1, 
          32}];

A21 = 0.50253532; A12 = 1.618147731;

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;

Total pressure in Pa

P = 101000;

Equilibrium relations

Table[
        y[i_, t_] := G1[i, t] PS1[i, t] x[i, t]/P, {i, 1, 
          32}];

EEq[i_] := 
      P ==  
        PS1[i, t] G1[i, t] x[i, t] + 
          PS2[i, t] G2[i, t] (1 - x[i, t])

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[EEq[i], {i, 1, 32}], 
          Table[x[i, 0] ==  0.5, {i, 1, 32}], 
          Table[T[i, 0] ==  350, {i, 1, 32}]], 
        Join[Table[x[i, t], {i, 1, 32}], Table[T[i, t], {i, 1, 32}]], {t, 0, 
          20000}]

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

Table[T[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[EEq[i], {i, 1, 32}], 
          Table[x[i, 0] ==  x[i, t] /. sol /. 
              t ->  20000, {i, 1, 32}], 
          Table[T[i, 0] ==  T[i, t] /. sol /. 
              t ->  20000, {i, 1, 32}]], 
        Join[Table[x[i, t], {i, 1, 32}], Table[T[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}]

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

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

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

Table[T[i, t] /. sol2 /. t ->  20000, {i, 1, 32}]