Application of the Riccati Equation to Solve Binary Distillation Problems​
​Authors: Housam Binous and Ahmed Bellagi

Acknowledgement :​
The Riccati equation obtained when solving a binary distillation problem can be found in
​Applied Mathematics and Modeling for Chemical Engineers by Rice and Do, Wiley, 1995.
We solve a binary distillation problem using the solution of a Riccati equation . We find that four plates are needed to achieve the desired product purity . We also compute the number of stages with the classical approach by stepping off stages using the equilibrium curve and the operating line . With this second approach, we find four plates too, in agreement with the Riccati equation method .
Solution of the Riccati equation
In[]:=
sol=RSolve[x[n+1]x[n]+Adx[n+1]+Bdx[n]+Cd0,x[n],n]//FullSimplify
Out[]=
x[n]-
n
4
Bd
n
-
1
Ad-Bd+
2
(Ad+Bd)
-4Cd
2
(Ad+Bd)
-4Cd
+
n
4
n
-
1
Ad-Bd+
2
(Ad+Bd)
-4Cd
(-Bd(Ad+Bd)+2Cd)+Bd
2
(Ad+Bd)
-4Cd
n
Ad-Bd+
2
(Ad+Bd)
-4Cd
AdBd-Cd

1
-(Bd(Ad+Bd)-2Cd)
n
Ad-Bd+
2
(Ad+Bd)
-4Cd
AdBd-Cd

1

n
4
(-Ad+Bd)
n
-
1
Ad-Bd+
2
(Ad+Bd)
-4Cd
+
n
4
n
-
1
Ad-Bd+
2
(Ad+Bd)
-4Cd
2
(Ad+Bd)
-4Cd
-(Ad-Bd)
n
Ad-Bd+
2
(Ad+Bd)
-4Cd
AdBd-Cd

1
-
2
(Ad+Bd)
-4Cd
n
Ad-Bd+
2
(Ad+Bd)
-4Cd
AdBd-Cd

1

Expressions for Ad, Bd and Cd for the distillation problem
In[]:=
Ad=
x0(α-1)-α(R+1)
(α-1)R
Out[]=
x0(-1+α)-(1+R)α
R(-1+α)
In[]:=
Bd=
1
(α-1)
Out[]=
1
-1+α
In[]:=
Cd=
x0
R(α-1)
Out[]=
x0
R(-1+α)
Relative volatility is equal to 2.45
In[]:=
α=2.45;
Reflux ratio is equal to 1.5143
In[]:=
R=1.5143;
Distillate mole fraction is equal to 0.9
In[]:=
x0=0.9;
Values of the constants Cd, Bd and Ad
In[]:=
Cd
Out[]=
0.409886
In[]:=
Bd
Out[]=
0.689655
In[]:=
Ad
Out[]=
-2.21112
We need exactly n=4 plates to obtain a mole fraction equal to 0.485144.
In[]:=
FindRoot[{First[(x[n]/.sol/.C[1]C1)]0.4851144,First[(x[n]/.sol/.C[1]C1/.n0)]x0},{C1,1},{n,5.5}]
Out[]=
{C13.62668,n4.00009}
Equilibrium function
In[]:=
equilb=Function[With[{a=2.45},(a#/(1+(a-1)#))]]
Out[]=
With{a=2.45},
a#1
1+(a-1)#1
&
Operating line
In[]:=
OpRect=Function[With[{R=1.5143},R/(R+1)#+x0/(R+1)]]
Out[]=
With{R=1.5143},
R#1
R+1
+
x0
R+1
&
Stepping off stages, we find the same value as above (i.e. n=4 plates)
In[]:=
y[0]=x0;x[0]=x0;i=0;
In[]:=
While[x[i]>0.48,{Print[x[i]," ",i," ",y[i]],v=FindRoot[y[i]equilb[u],{u,0.5}][[1,2]],t=OpRect[v],If[v>0.48,{x[i+1]=v,y[i+1]=t}],​​i++}]
0.9 0 0.9
0.786026 1 0.831356
0.668007 2 0.760276
0.564171 3 0.697738
0.48512 4 0.650128