Dynamics of a Five Equilibrium Stages 
Extraction Cascade with Solvent Recycle
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

First case: Solvent tank is a perfectly mixed container.
λ is proportional to the equilibrium constant. h depends on tank volume

Plotting the five stages compositions ( (Subscript[X, 1]), (Subscript[X, 2]), (Subscript[X, 3]), (Subscript[X, 4]), (Subscript[X, 5])) versus time for λ=1 and h=50.

λ = 1;

(Subscript[X, 0])[t] = 1;

h=50;

sol1 = 
      NDSolve[{(Subscript[X, 6])[t] - (1 + λ) 
                (Subscript[X, 5])[t] + λ (Subscript[X, 4])[t] == D[(Subscript[X, 5])[t], {t, 1}], 
          (Subscript[X, 5])[t] - (1 + λ) (Subscript[X, 4])[t] + λ (Subscript[X, 3])[t] == 
            D[(Subscript[X, 4])[t], {t, 1}], 
(Subscript[X, 4])[
                t] - (1 + λ) 
                (Subscript[X, 3])[t] + λ (Subscript[X, 2])[t] == D[(Subscript[X, 3])[t], {t, 1}], 
          (Subscript[X, 3])[t] - (1 + λ) (Subscript[X, 2])[t] + λ (Subscript[X, 1])[t] == 
            D[(Subscript[X, 2])[t], {t, 1}], 
          (Subscript[X, 2])[t] - (1 + λ) (Subscript[X, 1])[t] + λ (Subscript[X, 0])[t] == 
            D[(Subscript[X, 1])[t], {t, 1}], 
h D[
                (Subscript[X, 6])[t], {t, 1}] == ((Subscript[X, 1])[t] - (Subscript[X, 6])[t]), 
          (Subscript[X, 1])[0] == 0, (Subscript[X, 2])[0] == 0, (Subscript[X, 3])[0] == 0, 
          (Subscript[X, 4])[0] == 0, (Subscript[X, 5])[0] == 0, 
          (Subscript[X, 6])[0] == 0}, {(Subscript[X, 1])[t], (Subscript[X, 2])[t], (Subscript[X, 3])[t], (Subscript[X, 4])[t], (Subscript[X, 5])[t], 
          (Subscript[X, 6])[t]}, {t, 0, 150}]

Plot[
      Evaluate [
        Table[{(Subscript[X, 1])[t], (Subscript[X, 2])[t], (Subscript[X, 3])[t], (Subscript[X, 4])[t], (Subscript[X, 5])[t]}] /. sol1], {t, 0, 
        150}, PlotRange -> {0, 1}, 
      PlotStyle -> Table[Hue[(2 i - 1)/7], {i, 1, 5}]]

Plotting the first stage composition ( (Subscript[X, 1])) versus time for λ=1 (solid line) and 0.8 (dashed line).  h is set equal to 50, 100, 500 and 10^10.  We obtain the case without recycle for h very large (i.e. 10^10). For h small, we reach solvent saturation ( (Subscript[X, 1]) = 1) after a short period of time.

sol[h_, λ_] := 
      NDSolve[{(Subscript[X, 6])[t] - (1 + λ) 
                (Subscript[X, 5])[t] + λ (Subscript[X, 4])[t] == D[(Subscript[X, 5])[t], {t, 1}], 
          (Subscript[X, 5])[t] - (1 + λ) (Subscript[X, 4])[t] + λ (Subscript[X, 3])[t] == 
            D[(Subscript[X, 4])[t], {t, 1}], 
(Subscript[X, 4])[
                t] - (1 + λ) 
                (Subscript[X, 3])[t] + λ (Subscript[X, 2])[t] == D[(Subscript[X, 3])[t], {t, 1}], 
          (Subscript[X, 3])[t] - (1 + λ) (Subscript[X, 2])[t] + λ (Subscript[X, 1])[t] == 
            D[(Subscript[X, 2])[t], {t, 1}], 
          (Subscript[X, 2])[t] - (1 + λ) (Subscript[X, 1])[t] + λ (Subscript[X, 0])[t] == 
            D[(Subscript[X, 1])[t], {t, 1}], 
h D[
                (Subscript[X, 6])[t], {t, 1}] == ((Subscript[X, 1])[t] - (Subscript[X, 6])[t]), 
          (Subscript[X, 1])[0] == 0, (Subscript[X, 2])[0] == 0, (Subscript[X, 3])[0] == 0, 
          (Subscript[X, 4])[0] == 0, (Subscript[X, 5])[0] == 0, 
          (Subscript[X, 6])[0] == 0}, {(Subscript[X, 1])[t], (Subscript[X, 2])[t], (Subscript[X, 3])[t], (Subscript[X, 4])[t], (Subscript[X, 5])[t], 
          (Subscript[X, 6])[t]}, {t, 0, 150}]

plt1 = 
        Plot[Evaluate [
            Table[{(Subscript[X, 1])[t] /. sol[50, 1], (Subscript[X, 1])[t] /. sol[100, 1], 
                (Subscript[X, 1])[t] /. sol[500, 1], (Subscript[X, 1])[t] /. sol[10^10, 1]}]], {t, 0, 
            150}, PlotRange -> {0.5, 1}, 
          PlotStyle -> Table[Hue[(2 i - 1)/7], {i, 1, 5}], 
          Frame -> True, DisplayFunction -> Identity];

plt2 = 
        Plot[Evaluate [
            Table[{(Subscript[X, 1])[t] /. sol[50, 0.8], (Subscript[X, 1])[t] /. sol[100, 0.8], 
                (Subscript[X, 1])[t] /. sol[500, 0.8], (Subscript[X, 1])[t] /. sol[10^10, 0.8]}]], {t, 
            0, 150}, PlotRange -> {0.5, 1}, 
          PlotStyle -> {{Dashing[{0.02}], Hue[1/7]}, {Dashing[{0.02}], 
                Hue[3/7]}, {Dashing[{0.02}], Hue[5/7]}, {Dashing[{0.02}], 
                Hue[7/7]}}, Frame -> True, 
          DisplayFunction -> Identity];

Show[plt1,plt2,DisplayFunction->$DisplayFunction]

Second case: Solvent tank is a plug flow container.
λ is proportional to the equilibrium constant. h depends on tank volume

Plotting the first stage composition ( (Subscript[X, 1])) versus time for λ=1 and h=50.

λ = 1;

(Subscript[X, 0])[t] = 1;

h=50;

(Subscript[X, 6])[t] = 0;

sol2 = 
      NDSolve[{(Subscript[X, 6])[t] - (1 + λ) 
                (Subscript[X, 5])[t] + λ (Subscript[X, 4])[t] == D[(Subscript[X, 5])[t], {t, 1}], 
          (Subscript[X, 5])[t] - (1 + λ) (Subscript[X, 4])[t] + λ (Subscript[X, 3])[t] == 
            D[(Subscript[X, 4])[t], {t, 1}], 
(Subscript[X, 4])[
                t] - (1 + λ) 
                (Subscript[X, 3])[t] + λ (Subscript[X, 2])[t] == D[(Subscript[X, 3])[t], {t, 1}], 
          (Subscript[X, 3])[t] - (1 + λ) (Subscript[X, 2])[t] + λ (Subscript[X, 1])[t] == 
            D[(Subscript[X, 2])[t], {t, 1}], 
          (Subscript[X, 2])[t] - (1 + λ) (Subscript[X, 1])[t] + λ (Subscript[X, 0])[t] == 
            D[(Subscript[X, 1])[t], {t, 1}], (Subscript[X, 1])[0] == 0, (Subscript[X, 2])[0] == 0, 
          (Subscript[X, 3])[0] == 0, (Subscript[X, 4])[0] == 0, 
          (Subscript[X, 5])[0] == 0}, {(Subscript[X, 1])[t], (Subscript[X, 2])[t], (Subscript[X, 3])[t], (Subscript[X, 4])[t], 
          (Subscript[X, 5])[t]}, {t, 0, 50}]

plt3 = 
        Plot[Evaluate [
            Table[{(Subscript[X, 1])[t], (Subscript[X, 2])[t], (Subscript[X, 3])[t], (Subscript[X, 4])[t], (Subscript[X, 5])[t]}] /. 
              sol2], {t, 0, 50}, PlotRange -> {0, 1}, 
          PlotStyle -> Table[Hue[(2 i - 1)/7], {i, 1, 5}], 
          Frame -> True, DisplayFunction -> Identity];

(Subscript[X, 6])[t] = Flatten[Evaluate[(Subscript[X, 1])[t]] /. sol2][[1]]

sol3 = 
      NDSolve[{((Subscript[X, 6])[t] /. t -> t - 50) - (1 + λ) 
                (Subscript[X, 5])[t] + λ (Subscript[X, 4])[t] == D[(Subscript[X, 5])[t], {t, 1}], 
          (Subscript[X, 5])[t] - (1 + λ) (Subscript[X, 4])[t] + λ (Subscript[X, 3])[t] == 
            D[(Subscript[X, 4])[t], {t, 1}], 
(Subscript[X, 4])[
                t] - (1 + λ) 
                (Subscript[X, 3])[t] + λ (Subscript[X, 2])[t] == D[(Subscript[X, 3])[t], {t, 1}], 
          (Subscript[X, 3])[t] - (1 + λ) (Subscript[X, 2])[t] + λ (Subscript[X, 1])[t] == 
            D[(Subscript[X, 2])[t], {t, 1}], 
          (Subscript[X, 2])[t] - (1 + λ) (Subscript[X, 1])[t] + λ (Subscript[X, 0])[t] == 
            D[(Subscript[X, 1])[t], {t, 1}], 
          (Subscript[X, 1])[50] == Evaluate[(Subscript[X, 1])[t] /. sol2 /. 
                  t -> 50][[1]], 
          (Subscript[X, 2])[50] == Evaluate[(Subscript[X, 2])[t] /. sol2 /. 
                  t -> 50][[1]], 
          (Subscript[X, 3])[50] == Evaluate[(Subscript[X, 3])[t] /. sol2 /. 
                  t -> 50][[1]], 
          (Subscript[X, 4])[50] == Evaluate[(Subscript[X, 4])[t] /. sol2 /. 
                  t -> 50][[1]], 
          (Subscript[X, 5])[50] == Evaluate[(Subscript[X, 5])[t] /. sol2 /. 
                  t -> 50][[1]]}, {(Subscript[X, 1])[t], (Subscript[X, 2])[t], (Subscript[X, 3])[t], 
          (Subscript[X, 4])[t], (Subscript[X, 5])[t]}, {t, 50, 100}]

plt4 = 
        Plot[Evaluate [
            Table[{(Subscript[X, 1])[t], (Subscript[X, 2])[t], (Subscript[X, 3])[t], (Subscript[X, 4])[t], (Subscript[X, 5])[t]}] /. 
              sol3], {t, 50, 100}, PlotRange -> {0, 1}, 
          PlotStyle -> Table[Hue[(2 i - 1)/7], {i, 1, 5}], 
          Frame -> True, DisplayFunction -> Identity];

(Subscript[X, 6])[t] = Flatten[Evaluate[(Subscript[X, 1])[t]] /. sol3][[1]]

sol4 = 
      NDSolve[{((Subscript[X, 6])[t] /. t -> t - 50) - (1 + λ) 
                (Subscript[X, 5])[t] + λ (Subscript[X, 4])[t] == D[(Subscript[X, 5])[t], {t, 1}], 
          (Subscript[X, 5])[t] - (1 + λ) (Subscript[X, 4])[t] + λ (Subscript[X, 3])[t] == 
            D[(Subscript[X, 4])[t], {t, 1}], 
(Subscript[X, 4])[
                t] - (1 + λ) 
                (Subscript[X, 3])[t] + λ (Subscript[X, 2])[t] == D[(Subscript[X, 3])[t], {t, 1}], 
          (Subscript[X, 3])[t] - (1 + λ) (Subscript[X, 2])[t] + λ (Subscript[X, 1])[t] == 
            D[(Subscript[X, 2])[t], {t, 1}], 
          (Subscript[X, 2])[t] - (1 + λ) (Subscript[X, 1])[t] + λ (Subscript[X, 0])[t] == 
            D[(Subscript[X, 1])[t], {t, 1}], 
          (Subscript[X, 1])[100] == Evaluate[(Subscript[X, 1])[t] /. sol3 /. 
                  t -> 100][[1]], 
          (Subscript[X, 2])[100] == Evaluate[(Subscript[X, 2])[t] /. sol3 /. 
                  t -> 100][[1]], 
          (Subscript[X, 3])[100] == Evaluate[(Subscript[X, 3])[t] /. sol3 /. 
                  t -> 100][[1]], 
          (Subscript[X, 4])[100] == Evaluate[(Subscript[X, 4])[t] /. sol3 /. 
                  t -> 100][[1]], 
          (Subscript[X, 5])[100] == Evaluate[(Subscript[X, 5])[t] /. sol3 /. 
                  t -> 100][[1]]}, {(Subscript[X, 1])[t], (Subscript[X, 2])[t], (Subscript[X, 3])[t], 
          (Subscript[X, 4])[t], (Subscript[X, 5])[t]}, {t, 100, 150}]

plt5 = 
        Plot[Evaluate [
            Table[{(Subscript[X, 1])[t], (Subscript[X, 2])[t], (Subscript[X, 3])[t], (Subscript[X, 4])[t], (Subscript[X, 5])[t]}] /. 
              sol4], {t, 100, 150}, PlotRange -> {0, 1}, 
          PlotStyle -> Table[Hue[(2 i - 1)/7], {i, 1, 5}], 
          Frame -> True, DisplayFunction -> Identity];

Show[plt3,plt4,plt5,DisplayFunction->$DisplayFunction]