Ponchon Savarit Graphical 
Equilibrium Stage Method
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Suppression of warning messages

Off[FindRoot::lstol];

Off[General::spell1];

Off[InterpolatingFunction::dmval]

Vapor Pressures, Normal boiling points and total pressure 
Component 1= Hexane
Component 2= Octane

A2 = 6.4141; B2 = 5947.491; C2 = 360.26;
A1 = 6.039243; B1 = 5085.758; C1 = 382.794;

PS1 = 
        440 Exp[
            A1 - B1/(C1 + 1.8  T + 32)];
PS2 = 362.1 Exp[A2 - B2/(C2 + 1.8  T + 32)];

FindRoot[PS1 ==  14.696,  {T, 100}]

FindRoot[PS2 ==  14.696,  {T, 100}]

P = 14.696;

Activity coefficients using the Van Larr model

G1 = 
      Exp[A12 (A21 (1 - x)/(A12 x + A21 (1 - x)))^2] /. 
        x -> x[i]
G2 = 
      Exp[A21 (A12 x/(A12 x + A21 (1 - x)))^2] /. 
        x -> x[i]

A21 = -0.6731; A12 = -0.20100;

Vapor-liquid Equilibrium curve

h = 0.1; i = 0;

Do[{x[i] = h i, 
        T[i] = FindRoot[
              P ==  PS1 G1 x[i] + PS2 G2 (1 - x[i]), {T, 
                80}][[1, 2]], 
        y[i] =  PS1 G1 x[i]/P /. T -> T[i]}, {i, 0, 10}]

tbl1 = Table[{x[i], y[i]}, {i, 0, 10}]

plt1 = 
      ListPlot[tbl1, PlotJoined -> True, 
        PlotStyle -> RGBColor[1, 0, 0], 
        Epilog -> {RGBColor[0, 1, 0], 
            Line[{{0, 0}, {1, 1}}]}]

Saturated liquid and saturated vapor Enthalpy. h-x and H-y curves.

cl2[T_] := 2.85 10^2;
cl1[T_] := 2.29 10^2;

Do[hl[i] = 
        0.23901 (x[i] cl1[T[i] + 273.15] + (1 - x[i]) 
                cl2[T[i] + 273.15] ) T[i] , {i, 0, 10}]

tbl2 = Table[{x[i], hl[i]}, {i, 0, 10}]

G1 = 
      Exp[A12 (A21 (1 - x)/(A12 x + A21 (1 - x)))^2] /. 
        x -> xh
G2 = 
      Exp[A21 (A12 x/(A12 x + A21 (1 - x)))^2] /. 
        x -> xh

Do[{yh[i] = h i, 
        Th[i] = FindRoot[{P ==  PS1 G1 xh + PS2 G2 (1 - xh), 
                yh[i] ==   PS1 G1 xh/P}, {T, 100}, {xh, 
                0.2}][[1, 2]]}, {i, 0, 10}]

tbl3 = Table[{yh[i], Th[i]}, {i, 0, 10}]

CG1[T_] := -4.413 + 5.280 10^-1 T - 3.119 10^-4 T^2 + 
          6.494 10^-8 T^3;
CG2[T_] := -6.096 + 7.712 10^-1 T - 
          4.195  10^-4 T^2 + 8.855 10^-8 T^3;

dh1 = 
      2.87 10^4 ((507.4 - 273.15)/(507.4 - 68.9 - 
                  273.15))^0.38
dh2 = 
      3.44 10^4 ((568.8 - 273.15)/(568.8 - 125.6 - 
                  273.15))^0.38

Do[hv[i] = 
        0.23901 (yh[i] dh1 + (1 - yh[i]) 
                dh2 + (yh[i] CG1[Th[i] + 273.15] + (1 - yh[i]) 
                      CG2[Th[i] + 273.15] ) Th[i]) , {i, 0, 
        10}]

tbl4 = Table[{yh[i], hv[i]}, {i, 0, 10}]

plt2 = 
      ListPlot[tbl2, PlotJoined -> True, 
        PlotStyle -> RGBColor[1, 0, 0]]

plt3 = 
      ListPlot[tbl4, PlotJoined -> True, 
        PlotStyle -> RGBColor[0, 0, 1]]

Show[plt3, plt2, PlotRange -> {0, 18000}, 
      AspectRatio -> 1.5]

h-x and H-y curves with tie lines

vapor=Interpolation[tbl4,InterpolationOrder -> 1]

liquid=Interpolation[tbl2,InterpolationOrder -> 1]

tbl5 = 
      Table[Line[{{x[i], liquid[x[i]]}, {y[i], vapor[y[i]]}}], {i, 1, 
          9}]

tbl5bis = Table[{x[i], vapor[y[i]]}, {i, 0, 10}]

plt3bis = 
      ListPlot[tbl5bis, PlotStyle -> RGBColor[1, 0, 1], 
        PlotJoined -> True]

plt4 = 
      Show[plt3, plt2, plt3bis, Graphics[{RGBColor[0, 1, 0], tbl5}], 
        PlotRange -> {0, 18000}, AspectRatio -> 1.5, 
        Epilog -> {Line[{{x[2], liquid[x[2]]}, {x[2], 
                  vapor[y[2]]}}], 
Line[{{x[2], 
                  vapor[y[2]]}, {y[2], vapor[y[2]]}}]}]

Feed, distillate and Bottom compositions are 0.5, 0.95 and 0.05, respectively. Feed quality is q = 0.41. Feed enthalpy is hf.

xd[1] = 0.95; xb[1] = 0.05;

xf[1] = 0.5;

hf = Solve[
          0.41 ==  (vapor[xf[1]] - x)/(vapor[xf[1]] - 
                  liquid[xf[1]]), x][[1, 1, 2]]

Plotting Feed Point on Enthalpy composition diagram

Show[plt4, Graphics[{PointSize[0.03], Point[{xf[1], hf}]}]]

Interpolation to get tie lines

tbl6 = Table[{liquid[x[i]], vapor[y[i]]}, {i, 0, 10}]

tie1=Interpolation[tbl6,InterpolationOrder -> 1]

tbl7 = Table[{vapor[y[i]], liquid[x[i]]}, {i, 0, 10}]

tie2=Interpolation[tbl7,InterpolationOrder -> 1]

Computing the number of equilibrium stages at total reflux. The condenser is a partial condenser. 4 stages are required including the partial condenser and reboiler .

hd = vapor[xd[1]]

tie2[hd]

xe[1] = FindRoot[liquid[x] ==  tie2[hd], {x, 0.5}][[1, 
          2]]

ye[1] = xe[1]

xe[2] = FindRoot[
          liquid[x] ==  tie2[vapor[ye[1]]], {x, 0.5}][[1, 
          2]]

ye[2] = xe[2]

xe[3] = FindRoot[
          liquid[x] ==  tie2[vapor[ye[2]]], {x, 0.2}][[1, 
          2]]

ye[3] = xe[3]

xe[4] = FindRoot[
          liquid[x] ==  tie2[vapor[ye[3]]], {x, 0.2}][[1, 
          2]]

ye[0] = xd[1]

tbl8 = 
      Table[Line[{{ye[i], vapor[ye[i]]}, {xe[i + 1], 
              liquid[xe[i + 1]]}}], {i, 0, 3}]

tbl9 = 
      Table[Line[{{xe[i], liquid[xe[i]]}, {ye[i], vapor[ye[i]]}}], {i, 1, 
          3}]

plt4 = 
      Show[plt3, plt2, 
        Graphics[{RGBColor[0, 1, 0], tbl8, RGBColor[1, 0, 1], tbl9}], 
        PlotRange -> {0, 18000}, AspectRatio -> 1.5]

Computing difference points and L/D minimum.
We find L/D minimum = 1.03244.

FindRoot[{a (x1 - xf[1]) + hf ==  vapor[x1], 
        a (x2 - xf[1]) + hf ==  liquid[x2], 
        tie2[vapor[x1]] ==  liquid[x2]}, {a, 10000}, {x1, 0.5}, {x2, 
        0.3}]

line1[x_] := 18505.146261788144 (x - xf[1]) + hf

inters[1] = line1[0.95]

inters[2] = line1[0.05]

plt4 = 
      Show[plt3, plt2, 
        Graphics[{RGBColor[0, 1, 0], 
            Line[{{xd[1], inters[1]}, {xb[1], inters[2]}}], 
            RGBColor[1, 0, 1], PointSize[0.03], Point[{xd[1], inters[1]}], 
            Point[{xb[1], inters[2]}], Point[{xf[1], hf}]}], 
        PlotRange -> All, AspectRatio -> 1.5]

xe[1]

liquid[xe[1]]

xVn1 = FindRoot[(x1 - 
                    xd[1]) (liquid[xe[1]] - inters[1])/(xe[1] - 
                      xd [1]) + inters[1] ==  vapor[x1], {x1, 
            0.7}][[1, 2]]

LsurD = 
      Sqrt[(xVn1 - xd[1])^2 + (vapor[xVn1] - inters[1])^2]/
        Sqrt[(xVn1 - xe[1])^2 + (vapor[xVn1] - 
                  liquid[xe[1]])^2]

Computing difference Points when L/D real is 1.5 times greater than L/D minimum.

xe[1]

sol = 
      FindRoot[{a (x1 - xe[1]) + liquid[xe[1]] ==  
            vapor[x1], 
a (xd[1] - xe[1]) + 
              liquid[xe[1]] ==  
            intp, 
Sqrt[(x1 - xd[1])^2 + (vapor[x1] - 
                        intp)^2]/
              Sqrt[(x1 - xe[1])^2 + (vapor[x1] - 
                        liquid[xe[1]])^2] ==  1.5 LsurD}, {a, 
          80000}, {x1, 0.9}, {intp, 18000}]

ye[1] = sol[[2, 2]]

intp[1] = sol[[3, 2]]

intp[2] = (xb[1] - xf[1]) (intp[1] - hf)/(xd[1] - xf[1]) + 
        hf

plt4 = 
      Show[plt3, plt2, 
        Graphics[{RGBColor[0, 1, 0], 
            Line[{{xd[1], inters[1]}, {xb[1], inters[2]}}], 
            Line[{{xd[1], intp[1]}, {xb[1], intp[2]}}], RGBColor[1, 0, 1], 
            PointSize[0.03], Point[{xd[1], inters[1]}], 
            Point[{xb[1], inters[2]}], Point[{xd[1], intp[1]}], 
            Point[{xb[1], intp[2]}], Point[{xf[1], hf}]}], 
        PlotRange -> All, AspectRatio -> 1.5]

Stepping off stages in the rectifying section. We change operating line when the equilibrium tie line crosses the line joining the difference points. The optimum feed stage is stage 3 from the top including the partial condenser. The condenser is a partial condenser.

xe[2] = FindRoot[
          liquid[x] ==  tie2[vapor[ye[1]]], {x, 0.5}][[1, 
          2]]

sol = 
      FindRoot[{a (x1 - xe[2]) + liquid[xe[2]] ==  
            vapor[x1], 
a (xd[1] - xe[2]) + 
              liquid[xe[2]] ==  intp[1]}, {a, 100000}, {x1, 
          0.5}]

ye[2] = sol[[2, 2]]

xe[3] = FindRoot[
          liquid[x] ==  tie2[vapor[ye[2]]], {x, 0.4}][[1, 
          2]]

tbl10 = 
      Table[Line[{{ye[i], vapor[ye[i]]}, {xe[i + 1], 
              liquid[xe[i + 1]]}}], {i, 0, 2}]

\(\(\(tbl11 = 
      Table[Line[{{xe[i], liquid[xe[i]]}, {xd[1], intp[1]}}], {i, 1, 
          2}]\)\(
\)
    \)\)

plt4 = 
      Show[plt3, plt2, 
        Graphics[{RGBColor[1, 0, 1], 
            Line[{{xd[1], intp[1]}, {xb[1], intp[2]}}], Hue[0.35], tbl10, 
            Hue[0.88], tbl11, RGBColor[1, 0, 1], PointSize[0.03], 
            Point[{xd[1], intp[1]}], Point[{xb[1], intp[2]}], 
            Point[{xf[1], hf}]}], PlotRange -> All, 
        AspectRatio -> 1.5]

Stepping off stages in the stripping section. We change operating line when the equilibrium tie line crosses the line joining the difference points. The optimum feed stage is stage 3. The condenser is a partial condenser. A total of 7 stages is required including the partial condenser and reboiler.

sol = 
      FindRoot[{a (x1 - xe[3]) + liquid[xe[3]] ==  
            vapor[x1], 
a (xb[1] - xe[3]) + 
              liquid[xe[3]] ==  intp[2]}, {a, 100000}, {x1, 
          0.5}]

ye[3] = sol[[2, 2]]

xe[4] = FindRoot[
          liquid[x] ==  tie2[vapor[ye[3]]], {x, 0.4}][[1, 
          2]]

sol = 
      FindRoot[{a (x1 - xe[4]) + liquid[xe[4]] ==  
            vapor[x1], 
a (xb[1] - xe[4]) + 
              liquid[xe[4]] ==  intp[2]}, {a, 100000}, {x1, 
          0.5}]

ye[4] = sol[[2, 2]]

xe[5] = FindRoot[
          liquid[x] ==  tie2[vapor[ye[4]]], {x, 0.4}][[1, 
          2]]

sol = 
      FindRoot[{a (x1 - xe[5]) + liquid[xe[5]] ==  
            vapor[x1], 
a (xb[1] - xe[5]) + 
              liquid[xe[5]] ==  intp[2]}, {a, 100000}, {x1, 
          0.5}]

ye[5] = sol[[2, 2]]

xe[6] = FindRoot[
          liquid[x] ==  tie2[vapor[ye[5]]], {x, 0.4}][[1, 
          2]]

sol = 
      FindRoot[{a (x1 - xe[6]) + liquid[xe[6]] ==  
            vapor[x1], 
a (xb[1] - xe[6]) + 
              liquid[xe[6]] ==  intp[2]}, {a, 200000}, {x1, 
          0.1}]

ye[6] = sol[[2, 2]]

xe[7] = FindRoot[
          liquid[x] ==  tie2[vapor[ye[6]]], {x, 0.1}][[1, 
          2]]

tbl12 = 
      Table[Line[{{ye[i], vapor[ye[i]]}, {xe[i + 1], 
              liquid[xe[i + 1]]}}], {i, 3, 6}]

tbl13 = 
      Table[Line[{{ye[i], vapor[ye[i]]}, {xb[1], intp[2]}}], {i, 3, 
          6}]

plt4 = 
      Show[plt3, plt2, 
        Graphics[{RGBColor[0, 0, 0], 
            Line[{{xd[1], intp[1]}, {xb[1], intp[2]}}], Hue[0.35], tbl10, 
            tbl12, Hue[0.88], tbl11, Hue[0.58], tbl13, RGBColor[0, 0, 0], 
            PointSize[0.03], Point[{xd[1], intp[1]}], 
            Point[{xb[1], intp[2]}], Point[{xf[1], hf}]}], 
        PlotRange -> All, AspectRatio -> 1.5]

Plotting the McCabe and Thiele Diagram. Note that the operating lines are curved since we do not assume Equimolar Overflow. The feed line intersects with the operating lines approximately at x=0.345.

xe[0] = xd[1]

tbl14 = 
      Table[Line[{{xe[i - 1], ye[i - 1]}, {xe[i], ye[i - 1]}}], {i, 1, 
          7}]

tbl15 = 
      Table[Line[{{xe[i], ye[i - 1]}, {xe[i], ye[i]}}], {i, 1, 6}]

plt5 = Show[plt1, Graphics[{tbl14, tbl15}]]

tbl16 = Table[{xe[i], ye[i]}, {i, 0, 2}]

tbl17 = Table[{xe[i], ye[i]}, {i, 3, 6}]

op1=Interpolation[tbl16,InterpolationOrder -> 1]

op2=Interpolation[tbl17,InterpolationOrder -> 1]

Xint = FindRoot[op1[x] ==  op2[x], {x, 0.5}][[1, 
          2]]

Yint = op1[Xint]

plt6 = 
      ListPlot[Append[tbl16, {Xint, Yint}], PlotJoined -> True, 
        PlotStyle -> RGBColor[0, 0, 1]]

plt7 = 
      ListPlot[Sort[Append[tbl17, {Xint, Yint}]], PlotJoined -> True, 
        PlotStyle -> RGBColor[0, 0, 1]]

q = 1 - 0.59

feed[x_] := q/(q - 1) x - xf[1]/(q - 1)

plt8 = 
      Plot[feed[x], {x, 0, 0.5}, PlotStyle -> Hue[0.87]]

Show[plt5, plt6, plt7, plt8]

FindRoot[feed[x] == op2[x], {x, 0.3}]

FindRoot[feed[x] == op1[x], {x, 0.3}]

Computing the operating lines if we assume Equimolar Overflow. Comparing the operating lines with (pale blue) and without (dark blue) the Equimolar Overflow assumption.

DsurV = 1/(1 + LsurD)

OPRCMO[x_] := LsurD/(LsurD + 1) x + DsurV xd[1]

IntCMO[
        1] = FindRoot[feed[x] ==  OPRCMO[x], {x, 0.5}][[1, 
          2]]

IntCMO[2] = feed[IntCMO[1]]

Show[plt5, plt6, plt7, plt8, 
      Graphics[{Hue[0.55], 
          Line[{{xd[1], xd[1]}, {IntCMO[1], 
                IntCMO[2]}}], 
Line[{{xb[1], 
                xb[1]}, {IntCMO[1], IntCMO[2]}}]}]]