Application of the Riccati Equation 
to Solve Binary Distillation Problems
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

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.

Solution of the Riccati equation

sol = 
      RSolve[x[n + 1] x[n] + Ad x[n + 1] + Bd x[n] + Cd == 0, x[n], 
          n] // FullSimplify

Expressions for Ad, Bd and Cd for the distillation problem

Ad = ((x0 (α - 1) - α (R + 1))/((α - 1) R))

Bd = (1/((α - 1)))

Cd = (x0/(R (α - 1)))

Relative volatility is equal to 2.45

α = 2.45;

Reflux ratio is equal to 1.5143

R = 1.5143;

Distillate mole fraction is equal to 0.9

x0 = 0.9;

Values of the constants Cd, Bd and Ad

Cd

Bd

Ad

We need exactly n=4 plates to obtain a mole fraction equal to 0.485144.

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}]

Equilibrium function

equilb = 
      Function[With[{a = 
              2.45}, (a #/(1 + (a - 1) #))]]

Operating line

OpRect = 
      Function[With[{R = 1.5143}, 
          R/(R + 1) # + x0/(R + 1)]]

Stepping off stages, we find the same value as above (i.e. n=4 plates)

y[0] = x0; x[0] = x0; i = 0;

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++}]