Liquid chromatography with 
non-linear Langmuir adsorption equilibrium and transfer rate model    
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Mass Transfer coefficients of components A and B

keffA=1.88 10^-4;

keffB=1.94 10^-4;

Peclet number

Pe=0.2/ϵ+0.011/ϵ (Ren/ϵ)^0.48

Reynold number

Ren=ρ L (4/(ϵ Pi Dc^2)) Dp/η

Diffusivity

Di=L Dp 4/(Pe ϵ Dc^2 Pi)

Packing particle diameter=Dp
Column diameter=Dc
Void fraction of the bed=ϵ
Density of liquid phase=ρ
Viscosity of liquid phase=η
Volumetric flow rate=L

Dp=7.5 10^-3;Dc=4;L=0.4;ρ=0.74;η=2.2 10^-3;ϵ=0.4;

Injection of components A and B

(Subscript[CA, 0])[
        t_] := (Subscript[CA, init]) (UnitStep[t] - UnitStep[t - 10])

(Subscript[CA, init]) = 0.7;

Plot[(Subscript[CA, 0])[t], {t, 0, 200}, 
      PlotStyle -> RGBColor[1, 0, 0]]

(Subscript[CB, 0])[
        t_] := (Subscript[CB, init]) (UnitStep[t] - UnitStep[t - 10])

(Subscript[CB, init]) = 0.7;

Plot[(Subscript[CB, 0])[t], {t, 0, 200}, 
      PlotStyle -> RGBColor[1, 0, 0]]

Equilibrium constants of components A and B

(Subscript[K, A]) = 3;

(Subscript[K, B]) = 1.8;

Cross sectional area of the column

Ac=Pi Dc Dc/4

Area to volume ratio of the solid phase spheres

ap=6/Dp

Column length

Lbed=25;

length and number of a differential element

dz=Lbed/Ne

Ne=200;

Nonlinear Langmuir adsorption coefficients for components A and B

(Subscript[b, A]) = 2;  (Subscript[b, B]) = 4;

Solid phase equilibrium concentrations of components A and B. We assume non-linear Langmuir adsorption isotherms. Expressions of the tansfer rates labelled QA and QB.

(Subscript[CSAeq, 1])[t] = 
        (Subscript[K, A]) (Subscript[CA, 1])[
              t]/(1 + (Subscript[b, A]) (Subscript[CA, 1])[t] + (Subscript[b, B]) (Subscript[CB, 1])[t]);

(Subscript[CSBeq, 1])[t] = 
        (Subscript[K, B]) (Subscript[CB, 1])[
              t]/(1 + (Subscript[b, A]) (Subscript[CA, 1])[t] + (Subscript[b, B]) (Subscript[CB, 1])[t]);

(Subscript[CSAeq, Ne])[t] = 
        (Subscript[K, A]) (Subscript[CA, Ne])[
              t]/(1 + (Subscript[b, A]) (Subscript[CA, Ne])[t] + (Subscript[b, B]) (Subscript[CB, Ne])[t]);

(Subscript[CSBeq, Ne])[t] = 
        (Subscript[K, B]) (Subscript[CB, Ne])[
              t]/(1 + (Subscript[b, A]) (Subscript[CA, Ne])[t] + (Subscript[b, B]) (Subscript[CB, Ne])[t]);

i=2;

For[j = 2, 
      j < 200, {(Subscript[QA, i])[t] = 
          keffA ((Subscript[CSAeq, i])[t] - 
                (Subscript[CSA, i])[t]) (1 - ϵ) ap/ϵ, 
        (Subscript[CSAeq, i])[t] = 
          (Subscript[K, A]) (Subscript[CA, i])[t]/(1 + (Subscript[b, A]) (Subscript[CA, i])[t] + (Subscript[b, B]) (Subscript[CB, i])[t]), 
        i = i + 1}; j++]

i=2;

For[j = 2, 
      j < 200, {(Subscript[QB, i])[t] = 
          keffB ((Subscript[CSBeq, i])[t] - 
                (Subscript[CSB, i])[t]) (1 - ϵ) ap/ϵ, 
        (Subscript[CSBeq, i])[t] = 
          (Subscript[K, B]) (Subscript[CB, i])[t]/(1 + (Subscript[b, A]) (Subscript[CA, i])[t] + (Subscript[b, B]) (Subscript[CB, i])[t]), 
        i = i + 1}; j++]

(Subscript[QA, 1])[t] = 
        keffA ((Subscript[CSAeq, 1])[t] - 
              (Subscript[CSA, 1])[t]) (1 - ϵ) ap/ϵ;

(Subscript[QB, 1])[t] = 
        keffB ((Subscript[CSBeq, 1])[t] - 
              (Subscript[CSB, 1])[t]) (1 - ϵ) ap/ϵ;

(Subscript[QA, Ne])[t] = 
        keffA ((Subscript[CSAeq, Ne])[t] - 
              (Subscript[CSA, Ne])[t]) (1 - ϵ) ap/ϵ;

(Subscript[QB, Ne])[t] = 
        keffB ((Subscript[CSBeq, Ne])[t] - 
              (Subscript[CSB, Ne])[t]) (1 - ϵ) ap/ϵ;

Components balances for solid and liquid phases

eqA[1] = 
        D[(Subscript[CA, 1])[t], {t, 1}] == 
          Di ((Subscript[CA, 2])[t] - (Subscript[CA, 1])[t])/dz^2 - 
            L ((Subscript[CA, 2])[t] + (Subscript[CA, 1])[t] - 
                    2  (Subscript[CA, 0])[t])/(2 dz ϵ Ac) - 
            (Subscript[QA, 1])[t];

eqB[1] = 
        D[(Subscript[CB, 1])[t], {t, 1}] == 
          Di ((Subscript[CB, 2])[t] - (Subscript[CB, 1])[t])/dz^2 - 
            L ((Subscript[CB, 2])[t] + (Subscript[CB, 1])[t] - 
                    2  (Subscript[CB, 0])[t])/(2 dz ϵ Ac) - 
            (Subscript[QB, 1])[t];

eqA[i_] := 
      D[(Subscript[CA, i])[t], {t, 1}] == 
        Di ((Subscript[CA, i + 1])[t] - 2  (Subscript[CA, i])[t] + (Subscript[CA, i - 1])[t])/dz^2 - 
          L ((Subscript[CA, i + 1])[t] - 
                  (Subscript[CA, i - 1])[t])/(2 dz ϵ Ac) - 
          (Subscript[QA, i])[t]

eqB[i_] := 
      D[(Subscript[CB, i])[t], {t, 1}] == 
        Di ((Subscript[CB, i + 1])[t] - 2  (Subscript[CB, i])[t] + (Subscript[CB, i - 1])[t])/dz^2 - 
          L ((Subscript[CB, i + 1])[t] - 
                  (Subscript[CB, i - 1])[t])/(2 dz ϵ Ac) - 
          (Subscript[QB, i])[t]

eqSA[1] = 
        D[(Subscript[CSA, 1])[t], {t, 1}] == 
          keffA ap ((Subscript[CSAeq, 1])[t] - (Subscript[CSA, 1])[t]);

eqSB[1] = 
        D[(Subscript[CSB, 1])[t], {t, 1}] == 
          keffB ap ((Subscript[CSBeq, 1])[t] - (Subscript[CSB, 1])[t]);

eqSA[i_] := 
      D[(Subscript[CSA, i])[t], {t, 1}] == 
        keffA ap ((Subscript[CSAeq, i])[t] - (Subscript[CSA, i])[t])

eqSB[i_] := 
      D[(Subscript[CSB, i])[t], {t, 1}] == 
        keffB ap ((Subscript[CSBeq, i])[t] - (Subscript[CSB, i])[t])

eqSA[Ne] = 
        D[(Subscript[CSA, Ne])[t], {t, 1}] == 
          keffA ap ((Subscript[CSAeq, Ne])[t] - (Subscript[CSA, Ne])[t]);

eqSB[Ne] = 
        D[(Subscript[CSB, Ne])[t], {t, 1}] == 
          keffB ap ((Subscript[CSBeq, Ne])[t] - (Subscript[CSB, Ne])[t]);

eqA[Ne] = 
        D[(Subscript[CA, Ne])[t], {t, 1}] == 
          Di ((Subscript[CA, Ne - 1])[t] - (Subscript[CA, Ne])[t])/dz^2 - 
            L ((Subscript[CA, Ne])[t] - (Subscript[CA, Ne - 1])[t])/( 
                  dz ϵ Ac) - (Subscript[QA, Ne])[t];

eqB[Ne] = 
        D[(Subscript[CB, Ne])[t], {t, 1}] == 
          Di ((Subscript[CB, Ne - 1])[t] - (Subscript[CB, Ne])[t])/dz^2 - 
            L ((Subscript[CB, Ne])[t] - (Subscript[CB, Ne - 1])[t])/( 
                  dz ϵ Ac) - (Subscript[QB, Ne])[t];

Solving the set of differential equations

Unknown = 
        Join[Table[(Subscript[CSA, i]), {i, 1, 200}], Table[(Subscript[CSB, i]), {i, 1, 200}], 
          Table[(Subscript[CA, i]), {i, 1, 200}], Table[(Subscript[CB, i]), {i, 1, 200}]];

initialCond = 
        Join[Table[(Subscript[CA, i])[0] == 0, {i, 1, 200}], 
          Table[(Subscript[CB, i])[0] == 0, {i, 1, 200}], 
          Table[(Subscript[CSA, i])[0] == 0, {i, 1, 200}], 
          Table[(Subscript[CSB, i])[0] == 0, {i, 1, 200}]];

eqn=Join[Table[eqA[i],{i,1,200}],Table[eqB[i],{i,1,200}],Table[eqSA[i],{i,1,200}],Table[eqSB[i],{i,1,200}]];

sys=Flatten[Join[eqn,initialCond]];

sol=NDSolve[sys,Unknown,{t,0,5000}];

Plotting the profile of the concentration of components A and B at the exit of the column versus time. We get good separation of components A and B but both peaks present some skeweness because of the non-linear Langmuir equilibrium assumption.

plt1 = 
      Plot[Evaluate[(Subscript[CA, 200])[t]] /. sol, {t, 0, 3000}, PlotRange -> All, 
        PlotStyle -> RGBColor[1, 0, 1], 
        DisplayFunction -> Identity]

plt2 = 
      Plot[Evaluate[(Subscript[CB, 200])[t]] /. sol, {t, 0, 3000}, PlotRange -> All, 
        PlotStyle -> RGBColor[0, 1, 0], 
        DisplayFunction -> Identity]

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