Two components are separated by liquid chromatography using the two-phase rate model and either the stagewise model (finite differences) or orthogonal collocation. Linear equilibrium is assumed. The derivation of the governing equation are available in Ref. 1 where this problem was originally solved using Berkeley-Madonna and the stagewise method. Here, we solve the same problem using Mathematica and the orthogonal collocation. Finally, we compare the results obtained using both numerical techniques.
Mass Transfer coefficients of components A and B
In[]:=
keffA=1.88 10^-4;
In[]:=
keffB=1.94 10^-4;
Péclet number
In[]:=
Pe=0.2/ϵ+0.011/ϵ (Ren/ϵ)^0.48
Out[]=
0.2
ϵ
+
0.011
0.48

Ren
ϵ

ϵ
Reynolds number
In[]:=
Ren=ρ L (4/(ϵ Pi Dc^2)) Dp/η
Out[]=
4DpLρ
2
Dc
πϵη
Diffusivity
In[]:=
Di=L Dp 4/(Pe ϵ Dc^2 Pi)
Out[]=
4DpL
2
Dc
πϵ
0.2
ϵ
+
0.0123523
0.48
DpLρ
2
Dc
2
ϵ
η
ϵ
Nomenclature :​
Packing particle diameter=Dp
Column diameter=Dc
Volumetric flow rate=L
Density of liquid phase=ρ
Viscosity of liquid phase=η
Void fraction of the bed=ϵ
In[]:=
Dp=7.5 10^-3;
Dc=4;
L=0.4;
ρ=0.74;
η=2.2 10^-3;
ϵ=0.4;
Injection of components A and B
In[]:=
CA
0
[t_]:=
CA
init
(UnitStep[t]-UnitStep[t-10])
In[]:=
CA
init
=0.23;
In[]:=
Plot[
CA
0
[t],{t,0,200},PlotStyleRGBColor[1,0,0]]
Out[]=
In[]:=
CB
0
[t_]:=
CB
init
(UnitStep[t]-UnitStep[t-10])
In[]:=
CB
init
=0.23;
In[]:=
Plot[
CB
0
[t],{t,0,200},PlotStyleRGBColor[1,0,0]]
Out[]=
Equilibrium constants of components A and B
In[]:=
K
A
=3;
In[]:=
K
B
=1.8;
Cross sectional area of the column
In[]:=
Ac=Pi Dc Dc/4
Out[]=
4π
Area to volume ratio of the solid phase spheres
In[]:=
ap=6/Dp
Out[]=
800.
Column length
In[]:=
Lbed=25;
length and number of a differential element
In[]:=
dz=Lbed/Ne
Out[]=
25
Ne
In[]:=
Ne=200;

Ingham’s method based on finite differences

Solid phase equilibrium concentrations of components A and B. We assume linear adsorption isotherms. Expressions of the transfer rates labelled QA and QB.
In[]:=
CSAeq
1
[t]=
K
A
CA
1
[t];
In[]:=
CSBeq
1
[t]=
K
B
CB
1
[t];
In[]:=
CSAeq
Ne
[t]=
K
A
CA
Ne
[t];
In[]:=
CSBeq
Ne
[t]=
K
B
CB
Ne
[t];
In[]:=
For[i=2,i<200,{
QA
i
[t]=keffA(
CSAeq
i
[t]-
CSA
i
[t])(1-ϵ)ap/ϵ,
CSAeq
i
[t]=
K
A
CA
i
[t],i=i+1}];
In[]:=
For[i=2,i<200,{
QB
i
[t]=keffB(
CSBeq
i
[t]-
CSB
i
[t])(1-ϵ)ap/ϵ,
CSBeq
i
[t]=
K
B
CB
i
[t],i=i+1}];
In[]:=
QA
1
[t]=keffA(
CSAeq
1
[t]-
CSA
1
[t])(1-ϵ)ap/ϵ;
In[]:=
QB
1
[t]=keffB(
CSBeq
1
[t]-
CSB
1
[t])(1-ϵ)ap/ϵ;
In[]:=
QA
Ne
[t]=keffA(
CSAeq
Ne
[t]-
CSA
Ne
[t])(1-ϵ)ap/ϵ;
In[]:=
QB
Ne
[t]=keffB(
CSBeq
Ne
[t]-
CSB
Ne
[t])(1-ϵ)ap/ϵ;
Components balances for solid and liquid phases
In[]:=
eqA[1]=D[
CA
1
[t],{t,1}]Di(
CA
2
[t]-
CA
1
[t])/dz^2-L(
CA
2
[t]+
CA
1
[t]-2
CA
0
[t])/(2dzϵAc)-
QA
1
[t];
In[]:=
eqB[1]=D[
CB
1
[t],{t,1}]Di(
CB
2
[t]-
CB
1
[t])/dz^2-L(
CB
2
[t]+
CB
1
[t]-2
CB
0
[t])/(2dzϵAc)-
QB
1
[t];
In[]:=
eqA[i_]:=D[
CA
i
[t],{t,1}]Di(
CA
i+1
[t]-2
CA
i
[t]+
CA
i-1
[t])/dz^2-L(
CA
i+1
[t]-
CA
i-1
[t])/(2dzϵAc)-
QA
i
[t]
In[]:=
eqB[i_]:=D[
CB
i
[t],{t,1}]Di(
CB
i+1
[t]-2
CB
i
[t]+
CB
i-1
[t])/dz^2-L(
CB
i+1
[t]-
CB
i-1
[t])/(2dzϵAc)-
QB
i
[t]
In[]:=
eqSA[1]=D[
CSA
1
[t],{t,1}]keffAap(
CSAeq
1
[t]-
CSA
1
[t]);
In[]:=
eqSB[1]=D[
CSB
1
[t],{t,1}]keffBap(
CSBeq
1
[t]-
CSB
1
[t]);
In[]:=
eqSA[i_]:=D[
CSA
i
[t],{t,1}]keffAap(
CSAeq
i
[t]-
CSA
i
[t])
In[]:=
eqSB[i_]:=D[
CSB
i
[t],{t,1}]keffBap(
CSBeq
i
[t]-
CSB
i
[t])
In[]:=
eqSA[Ne]=D[
CSA
Ne
[t],{t,1}]keffAap(
CSAeq
Ne
[t]-
CSA
Ne
[t]);
In[]:=
eqSB[Ne]=D[
CSB
Ne
[t],{t,1}]keffBap(
CSBeq
Ne
[t]-
CSB
Ne
[t]);
In[]:=
eqA[Ne]=D[
CA
Ne
[t],{t,1}]Di(
CA
Ne-1
[t]-
CA
Ne
[t])/dz^2-L(
CA
Ne
[t]-
CA
Ne-1
[t])/(dzϵAc)-
QA
Ne
[t];
In[]:=
eqB[Ne]=D[
CB
Ne
[t],{t,1}]Di(
CB
Ne-1
[t]-
CB
Ne
[t])/dz^2-L(
CB
Ne
[t]-
CB
Ne-1
[t])/(dzϵAc)-
QB
Ne
[t];
Solving the set of differential equations
In[]:=
Unknown=Join[Table[
CSA
i
,{i,1,200}],Table[
CSB
i
,{i,1,200}],Table[
CA
i
,{i,1,200}],Table[
CB
i
,{i,1,200}]];
In[]:=
initialCond=Join[Table[
CA
i
[0]0,{i,1,200}],Table[
CB
i
[0]0,{i,1,200}],Table[
CSA
i
[0]0,{i,1,200}],Table[
CSB
i
[0]0,{i,1,200}]];
In[]:=
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}]];
In[]:=
sys=Flatten[Join[eqn,initialCond]];
In[]:=
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. Peaks are quite symmetric because of the linear equilibrium assumption. Peak two is wider because its elution time is larger.
In[]:=
plt1b=Plot[Evaluate[
CA
200
[t]]/.sol,{t,0,3000},PlotRangeAll,PlotStyle{RGBColor[0,0,0],Dashed}];
In[]:=
plt2b=Plot[Evaluate[
CB
200
[t]]/.sol,{t,0,3000},PlotRangeAll,PlotStyle{Dashed,RGBColor[0,0,0]}];
In[]:=
Show[plt1b,plt2b,ImageSize-> Medium,AspectRatio-> 1,Frame-> True,GridLines->Automatic,
FrameLabel->{Style["time",14],Style["liquid concentration of A and B",14]}]
Out[]=
Our method based on orthogonal collocation and the Gauss-Lobatto nodes
In[]:=
ClearAll["CA*","CB*","QA[*","QB[*","CSA*","CSB*","eqA*","eqB*","eqS*"]
Chebyshev Derivative matrices (41 collocation points)
In[]:=
Np=40;​​d[0,0]=(2Np^2+1)/6//N;​​d[Np,Np]=-(2Np^2+1)/6//N;​​yc[i_]:=Cos[iPi/Np]//N;​​d[i_,i_]:=-yc[i]/(2(1-yc[i]^2))//N;​​d[j_,k_]:=If[j≠k&&k≠0&&k≠Np&&j≠0&&j≠Np,(-1)^(j+k)/(yc[j]-yc[k])]//N;​​d[Np,0]=-1/2(-1)^Np;​​d[0,Np]=1/2(-1)^Np;​​d[j_,0]:=If[j≠Np&&j≠0,-1/2(-1)^j/(1-yc[j])];​​d[j_,Np]:=If[j≠Np&&j≠0,1/2(-1)^(j+Np)/(1+yc[j])];​​d[0,j_]:=If[j≠Np&&j≠0,2(-1)^j/(1-yc[j])];​​d[Np,j_]:=If[j≠Np&&j≠0,-2(-1)^(j+Np)/(1+yc[j])];​​Add=Table[d[i,j]//N,{i,0,Np},{j,0,Np}];Bdd=Add.Add;
First-order and Second-order derivatives of CA and CB
Solid phase equilibrium concentrations of components A and B. We assume linear adsorption isotherms. Expressions of the transfer rates labelled QA and QB.
Components balances for solid and liquid phases
Solving the set of differential equations
In[]:=
eqn=Join[Table[eqA[i],{i,1,Np}],Table[eqB[i],{i,1,Np}],Table[eqSA[i],{i,1,Np+1}],Table[eqSB[i],{i,1,Np+1}]];
In[]:=
sys=Flatten[Join[eqn,initialCond]];
In[]:=
sol=NDSolve[sys,Unknown,{t,0,5000},AccuracyGoal -> 10];
Plotting the profile of the concentration of components A and B at the exit of the column versus time. Peaks are quite symmetric because of the linear equilibrium assumption. Peak two is wider because its elution time is larger.
In[]:=
Show[plt1,AspectRatio-> 1,Frame-> True,GridLines->Automatic,
FrameLabel->{Style["time",14],Style["liquid concentration of A and B",14]}]
Comparison between our method based on orthogonal colocation (colored lines) and Ingham' s method based on finite differences (black dashed lines)
In[]:=
Show[plt1,plt1b,plt2b,AspectRatio-> 1,Frame-> True,GridLines->Automatic,
FrameLabel->{Style["time",14],Style["liquid concentration of A and B",14]}]

References

◼
  • [1] J. Ingham, I. J. Dunn, E. Heinzle, J. E. Prenosil, J.B. Snape, Chemical Engineering Dynamics, Weinheim: Wiley - VCH, 2000

    • CITE THIS NOTEBOOK

    Liquid Chromatography with Linear Adsorption Equilibrium and Transfer Rate Model​
    by Housam Binous and Ahmed Bellagi
    Wolfram Community, STAFF PICKS, March 19, 2024
    ​https://community.wolfram.com/groups/-/m/t/3143304