Four stage batch distillation column
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

We use a batch column with four theoretical plates to separate a binary mixture of water and methanol at 1 bar. Temperature and concentration profiles versus time are obtained for a reflux ratio equal to 10. In addition, the total moles in the still and distillate are plotted versus time. This demo shows how one can use NDSolve to treat a complex system of 13 differential algebraic equations very easily.

Reference: Example BUBBLE in page 559 by J. Ingham, I. J. Dunn, E. Heinzle and J. E. Prenosil, Chemical Engineering Dynamics, Wiley-VCH, Weinheim, Second Edition, 2000.

Message Suppression

Off[General::spell1]

Component balances for the methanol on the plates, the still and condenser

diffeq1=M D[x1[t],{t,1}]==L (x0[t]-x1[t])+V (y2[t]-y1[t])

diffeq2=M D[x2[t],{t,1}]==L (x1[t]-x2[t])+V (y3[t]-y2[t])

diffeq3=M D[x3[t],{t,1}]==L (x2[t]-x3[t])+V (y4[t]-y3[t])

diffeq4=M D[x4[t],{t,1}]==L (x3[t]-x4[t])+V (yB[t]-y4[t])

diffeq5=D[MB[t],{t,1}]==L -V

diffeq0=M0 D[x0[t],{t,1}]==V y1[t]-(L+Di) x0[t]

diffeq6=D[MB[t] xB[t],{t,1}]==L x4[t]-V yB[t]

diffeq7=D[DIST[t],{t,1}]==Di

Total mass balance equation and reflux ratio definition

Di=V/(R+1);

L=V-Di;

Initial conditions

init0=x0[0]==0;

init1=x1[0]==0;

init2=x2[0]==0;

init3=x3[0]==0;

init4=x4[0]==0;

initMB=MB[0]==2500;

initDIST=DIST[0]==0;

initXB=xB[0]==0.8;

init5=T1[0]==50;

init6=T2[0]==50;

init7=T3[0]==50;

init8=T4[0]==50;

init9=TB[0]==50;

Values of parameters used in the expressions to compute the activity coefficients and vapor pressures. We take the Van Laar equation to estimate the activity coefficients since the binary mixture: water-methanol presents deviation from ideality in the liquid phase.

D1=-4617.8;

C1=13.676;

D2=-5042.6;

C2=13.519;

A=0.85;

B=0.48;

Values of the reflux ratio, the vapor flow rate, the molar holdup of the plates and the reflux drum

R=10;

V=10;

M0=100;

M=20;

Total pressure is set equal to 1 bar

P=1;

Algebraic relations representing VL Equilibrium

equilb1=P==Exp[A*(1-x1[t]^2)/(1-x1[t]+A x1[t]/B)^2] Exp[C1+D1/(273.15+T1[t])] x1[t]+
Exp[A*x1[t]^2/(x1[t]+A (1-x1[t])/B)^2] Exp[C2+D2/(273.15+T1[t])] (1-x1[t]);

equilb2=P==Exp[A*(1-x2[t]^2)/(1-x2[t]+A x2[t]/B)^2] Exp[C1+D1/(273.15+T2[t])] x2[t]+
Exp[A*x2[t]^2/(x2[t]+A (1-x2[t])/B)^2] Exp[C2+D2/(273.15+T2[t])] (1-x2[t]);

equilb3=P==Exp[A*(1-x3[t]^2)/(1-x3[t]+A x3[t]/B)^2] Exp[C1+D1/(273.15+T3[t])] x3[t]+
Exp[A*x3[t]^2/(x3[t]+A (1-x3[t])/B)^2] Exp[C2+D2/(273.15+T3[t])] (1-x3[t]);

equilb4=P==Exp[A*(1-x4[t]^2)/(1-x4[t]+A x4[t]/B)^2] Exp[C1+D1/(273.15+T4[t])] x4[t]+
Exp[A*x4[t]^2/(x4[t]+A (1-x4[t])/B)^2] Exp[C2+D2/(273.15+T4[t])] (1-x4[t]);

equilbB=P==Exp[A*(1-xB[t]^2)/(1-xB[t]+A xB[t]/B)^2] Exp[C1+D1/(273.15+TB[t])] xB[t]+
Exp[A*xB[t]^2/(xB[t]+A (1-xB[t])/B)^2] Exp[C2+D2/(273.15+TB[t])] (1-xB[t]);

y1[t_]:=Exp[A*(1-x1[t]^2)/(1-x1[t]+A x1[t]/B)^2] Exp[C1+D1/(273.15+T1[t])] x1[t]/P;

y2[t_]:=Exp[A*(1-x2[t]^2)/(1-x2[t]+A x2[t]/B)^2] Exp[C1+D1/(273.15+T2[t])] x2[t]/P;

y3[t_]:=Exp[A*(1-x3[t]^2)/(1-x3[t]+A x3[t]/B)^2] Exp[C1+D1/(273.15+T3[t])] x3[t]/P;

y4[t_]:=Exp[A*(1-x4[t]^2)/(1-x4[t]+A x4[t]/B)^2] Exp[C1+D1/(273.15+T4[t])] x4[t]/P;

yB[t_]:=Exp[A*(1-xB[t]^2)/(1-xB[t]+A xB[t]/B)^2] Exp[C1+D1/(273.15+TB[t])] xB[t]/P;

The number of the differential algebraic equations is 13

Length[{diffeq1,diffeq2,diffeq3,diffeq4,diffeq5,diffeq0,diffeq6,diffeq7,equilb1,equilb2,equilb3,equilb4,equilb5}]

The number of the initial conditions is 13

Length[{init0,init1,init2,init3,init4,init5,init6,init7,init8,init9,initXB,initMB,initDIST}]

The number of the functions is 13

Length[{x0[t],x1[t],x2[t],x3[t],x4[t],xB[t],T1[t],T2[t],T3[t],T4[t],TB[t],MB[t],DIST[t]}]

Calling NDSolve to solve the system of 13 DAEs

sol=NDSolve[{diffeq1,diffeq2,diffeq3,diffeq4,diffeq5,diffeq0,diffeq6,diffeq7,equilb1,equilb2,equilb3,equilb4,equilbB,init0,init1,init2,init3,init4,init5,init6,init7,init8,init9,initXB,initMB,initDIST},{x0[t],x1[t],x2[t],x3[t],x4[t],xB[t],T1[t],T2[t],T3[t],T4[t],TB[t],MB[t],DIST[t]},{t,0,200}]

Plotting the total number of moles in the still

plt1=Plot[Evaluate[MB[t]/.sol],{t,0,200},PlotStyle->RGBColor[1,0,0],PlotRange->All]

Plotting the concentration profile versus time in the still, the plates and the condenser

plt2=Plot[Evaluate[xB[t]/.sol],{t,0,200},PlotStyle->RGBColor[0,1,0],PlotRange->All]

plt3=Plot[Evaluate[x0[t]/.sol],{t,0,200},PlotStyle->RGBColor[1,0,0],PlotRange->All]

plt4=Plot[Evaluate[x1[t]/.sol],{t,0,200},PlotStyle->RGBColor[0,0,1],PlotRange->All]

plt5=Plot[Evaluate[x2[t]/.sol],{t,0,200},PlotStyle->RGBColor[1,0,1],PlotRange->All]

plt6=Plot[Evaluate[x3[t]/.sol],{t,0,200},PlotStyle->RGBColor[0.75,1,0],PlotRange->All]

plt7=Plot[Evaluate[x4[t]/.sol],{t,0,200},PlotStyle->RGBColor[1,0,0.75],PlotRange->All]

Show[plt2,plt3,plt4,plt5,plt6,plt7]

Plotting the temperature profile versus time in the still and in the plates

plt8=Plot[Evaluate[T1[t]/.sol],{t,0,200},PlotStyle->RGBColor[1,0,0],PlotRange->All]

plt9=Plot[Evaluate[T2[t]/.sol],{t,0,200},PlotStyle->RGBColor[0,0,1],PlotRange->All]

plt10=Plot[Evaluate[T3[t]/.sol],{t,0,200},PlotStyle->RGBColor[1,0,1],PlotRange->All]

plt11=Plot[Evaluate[T4[t]/.sol],{t,0,200},PlotStyle->RGBColor[0.75,1,0],PlotRange->All]

plt12=Plot[Evaluate[TB[t]/.sol],{t,0,200},PlotStyle->RGBColor[0.75,0.75,0],PlotRange->All]

Show[plt8,plt9,plt10,plt11,plt12,PlotRange->{60,100}]

Plotting the total number of moles in the distillate

plt13=Plot[Evaluate[DIST[t]/.sol],{t,0,200},PlotStyle->RGBColor[0.75,1,0.5],PlotRange->All]