Consider a distillation column operating at 101.325 kPa and separating an equimolar mixture of benzene (
B
), toluene (
T
), and
p
-xylene (
X
). The column has 20 stages, a partial reboiler, and a total condenser. The feed location is stage 13 counting from the top. The feed, a saturated liquid, has a flow rate equal to 100 kmole/hr. The column has a liquid sidestream at stage 7 counting from the top. You can set the draw rate of this side-stream as well as the reflux and reboil ratios of the column. The We plots the composition profile in orange, green, and blue for benzene, toluene, and
p
-xylene, respectively. Plots of the temperature profile and liquid and vapor flow rates are also given. The variables
F
,
D
,
R
, and
S
stand for feed, distillate, residue, and sidestream flow rates, respectively. The column’s degree of freedom is 3 and we choose to fix: the reflux ratio, the reboil ratio and the sidestream draw rate. Finally, we compare with the results of PYTHON and provide the corresponding PYTHON coding.
In[]:=
Srate=15;RReflux=5;RReboil=5;​​​​(*ModifiedAntoineConstantsforBenzene-toluene-p-xylenereproducedfromAspen-HYSYS*)​​A1=169.65;B1=-10314.8;C1=0.;D1=-23.5895;E1=
2.09442
5
10
;F1=2.;A2=76.4511;B2=-6995.;C2=0.;D2=-9.1635;E2=
6.225
6
10
;F2=2.;A3=91.2505;B3=-8121.25;C3=0.;D3=-11.3188;E3=
7.30555
6
10
;F3=2.;​​​​PSat[1][i_]:=
E1
F1
T
+A1+D1Log[T]+
B1
C1+T

/.TT[i];PSat[2][i_]:=
E2
F2
T
+A2+D2Log[T]+
B2
C2+T

/.TT[i];PSat[3][i_]:=
E3
F3
T
+A3+D3Log[T]+
B3
C3+T

/.TT[i];
Constants for liquid and vapor enthalpies for Benzene - toluene - p-xylene reproduced from Aspen-HYSYS
a1=0.;b1=129.44;c1=-
8.475001
2
10
;d1=
2.159366
4
10
;e1=-
2.089645
14
10
;h1=e1
4
T
+d1
3
T
+c1
2
T
+b1T+a1/.TTF;hl1=e1
4
T
+d1
3
T
+c1
2
T
+b1T+a1-h1;a2=0.;b2=140.1399;c2=-
7.615016
2
10
;d2=
2.31667
4
10
;e2=-
2.769433
14
10
;h2=e2
4
T
+d2
3
T
+c2
2
T
+b2T+a2/.TTF;hl2=e2
4
T
+d2
3
T
+c2
2
T
+b2T+a2-h2;a3=0.;b3=-35.49996;c3=0.6435998;d3=-
8.66333
4
10
;e3=
6.064997
7
10
;h3=e3
4
T
+d3
3
T
+c3
2
T
+b3T+a3/.TTF;hl3=e3
4
T
+d3
3
T
+c3
2
T
+b3T+a3-h3;​​​​a1v=45880.06;b1v=114.7825;c1v=-0.2426732;d1v=
5.663124
4
10
;e1v=-
3.654795
7
10
;Hv1=e1v
4
T
+d1v
3
T
+c1v
2
T
+b1vT+a1v-h1;a2v=48522.49;b2v=132.693;c2v=-0.2388932;d2v=
5.907287
4
10
;e2v=-
3.913967
7
10
;Hv2=e2v
4
T
+d2v
3
T
+c2v
2
T
+b2vT+a2v-h2;a3v=11653.17;b3v=383.7983;c3v=-1.17436;d3v=
2.308164
3
10
;e3v=-
1.543837
6
10
;Hv3=e3v
4
T
+d3v
3
T
+c3v
2
T
+b3vT+a3v-h3;
Liquid\b-phase and vapor-phase enthalpies for mixtures of Benzene - toluene - p-xylene
H[i_]:=Hv1y[1][i]+Hv2y[2][i]+Hv3y[3][i]/.TT[i];​​h[i_]:=hl1x[1][i]+hl2x[2][i]+hl3x[3][i]/.TT[i];
Feed thermal quality and composition
z[i_]:=1/3;phi=0;
Total pressure
P=101.325;
Some calculation regarding the feed stream
solf=QuietFindRoot
3
∑
i=1

PSat[i][1]
P
-1z[i]

PSat[i][1]
P
-1phi+1
0,{T[1],350,420};​​​​TF=T[1]/.solf(*feedstreamtemeprature*)​​​​hF=(1-phi)(hl1xz[1]+hl2xz[2]+hl3xz[3])+phi(Hv1yz[1]+Hv2yz[2]+Hv3yz[3])/.TTF;(*feedstreamentahlpy*)​​​​hF​​​​(*vaporinequilibrium---itscomposition*)​​yz[1]=
PSat[1][1]xz[1]
P
/.solf;yz[2]=
PSat[2][1]xz[2]
P
/.solf;yz[3]=
PSat[3][1]xz[3]
P
/.solf;yz[4]=
PSat[4][1]xz[4]
P
/.solf;​​​​(*liquidinequilibrium---itscomposition*)​​xz[1]=
z[1]

PSat[1][1]
P
-1phi+1
/.solf;xz[2]=
z[2]

PSat[2][1]
P
-1phi+1
/.solf;xz[3]=
z[3]

PSat[3][1]
P
-1phi+1
/.solf;xz[4]=
z[4]

PSat[4][1]
P
-1phi+1
/.solf;​​​​Table[xz[i],{i,1,3}]​​​​Table[yz[i],{i,1,3}]
Out[]=
374.497
Out[]=
3.03165×
-12
10
Out[]=

1
3
,
1
3
,
1
3

Out[]=
{0.634953,0.254035,0.111012}
Summation rules
y[3][i_]:=-y[1][i]-y[2][i]+1;x[3][i_]:=-x[1][i]-x[2][i]+1;
Feed flow-rate
F=100;​​(*MESHEquations*)​​eq1=Di+L[0]V[1];​​eq2=Dih[0]+L[0]h[0]V[1]H[1]+Qc;​​​​eq3[i_]:=V[i]+L[i]+If[i7,S,0]If[i13,F,0]+V[i+1]+L[i-1];​​eq4[i_]:=Table[V[i]y[j][i]+L[i]x[j][i]+If[i7,Sx[j][i],0]If[i13,Fz[j],0]+V[i+1]y[j][i+1]+L[i-1]x[j][i-1],{j,1,2}];​​eq5[i_]:=Table[Py[j][i]PSat[j][i]x[j][i],{j,1,3}];​​eq6[i_]:=V[i]H[i]+L[i]h[i]+If[i7,Sh[i],0]If[i13,FhF,0]+V[i+1]H[i+1]+L[i-1]h[i-1];​​​​eq15=V[21]+BL[20];​​eq16=Table[V[21]y[i][21]+Bx[i][21]L[20]x[i][20],{i,1,2}];​​eq17=Table[Py[i][21]PSat[i][21]x[i][21],{i,1,3}];​​eq18=V[21]H[21]+Bh[21]L[20]h[20]+Qr;​​​​eq19=V[21]/BRReboil;​​eq20=L[0]/DiRReflux;​​eq21=SSrate;​​​​eq22=Table[x[i][0]y[i][1],{i,1,2}];​​​​eq25=Table[Py[i][0]PSat[i][0]x[i][0],{i,1,3}];
Solving MESH equations with FindRoot
Plotting results
PLT=ListPlot[{Table[{i,x[1][i]}/.sol1e,{i,0,21}],Table[{i,x[2][i]}/.sol1e,{i,0,21}],Table[{i,x[3][i]}/.sol1e,{i,0,21}]},PlotMarkers{"◆",16},AxesOrigin{0,0},JoinedTrue,FrameLabel{Style["plate number ",16],Style["steady-state composition ",16]},FrameTrue,GridLinesAutomatic,PlotStyle{{Thick,Orange},{Thick,RGBColor[0,1,0]},{Thick,Blue}},ImageSize{350,350},AspectRatio1,PlotRange{-0.05,1.05}]​​
Out[]=
In[]:=
ListPlot[{Table[{i,-273.15+T[i]}/.sol1e,{i,0,21}]},PlotMarkers{"◆",16},AxesOrigin{0,0},JoinedTrue,FrameLabel{Style["plate number ",16],Style["steady-state temperature ",16]},FrameTrue,GridLinesAutomatic,PlotStyle{{Thick,Orange},{Thick,RGBColor[0,1,0]},{Thick,Blue}},ImageSize{350,350},AspectRatio1,PlotRange{70,150}]​​
Out[]=
Sketch of the distillation column
Python code is attached at the end.

References

[1] E. J. Henley and J. D. Seader, Equilibrium-Stage Separation Operations in Chemical Engineering, New York: Wiley, 1981.

CITE THIS NOTEBOOK

Multicomponent distillation column with a sidestream​
by Housam Binous and Ahmed Bellagi
Wolfram Community, STAFF PICKS, October 20, 2025
​https://community.wolfram.com/groups/-/m/t/3562838