Two Feed Extraction Column
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Liquid Liquid Equilibrium of ternary system Isopropyl ether-water-isopropyl alcohol.
Component 1 = water 
Component 2= Acetic Acid
Component 3=Isopropyl ether

tbl1 = {{98.1, 0.69}, {97.1, 1.41}, {95.5, 2.89}, {91.7, 
            6.42}, {84.4, 13.30}, {71.1, 25.50}, {58.9, 36.70}, {45.1, 
            44.30}, {37.1, 46.40}};

tbl2 = {{0.5, 0.18}, {0.7, 0.37}, {0.8, 0.79}, {1.0, 1.93}, {1.9, 
            4.82}, {3.9, 11.40}, {6.9, 21.60}, {10.8, 31.10}, {15.1, 
            36.20}};

plt1 = 
      ListPlot[tbl1, PlotJoined -> True, 
        PlotStyle -> RGBColor[1, 0, 0]]

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

Show[plt1, plt2, PlotRange -> {{0, 100}, {0, 100}}, 
      Epilog -> Line[{{100, 0}, {0, 100}}]]

rafnt=Interpolation[tbl1,InterpolationOrder -> 1]

rafnt[90]

Extrc=Interpolation[tbl2,InterpolationOrder -> 1]

plt3 = 
      Plot[rafnt[x], {x, 37.1, 98.1}, PlotRange -> {0, 100}, 
        PlotStyle -> RGBColor[1, 0, 0]]

plt4 = 
      Plot[Extrc[x], {x, 0.5, 15.1}, PlotRange -> {0, 100}, 
        PlotStyle -> RGBColor[0, 0, 1]]

Show[plt3, plt4, PlotRange -> {{0, 100}, {0, 100}}, 
      Epilog -> Line[{{100, 0}, {0, 100}}]]

Using interpolation to get tie lines

tbl3 = {{0.5, 98.1}, {0.7, 97.1}, {0.8, 95.5}, {1.0, 91.7}, {1.9, 
            84.4}, {3.9, 71.1}, {6.9, 58.9}, {10.8, 45.1}, {15.1, 
            37.1}};

tie1=Interpolation[tbl3,InterpolationOrder -> 1]

tbl4 = {{98.1, 0.5}, {97.1, 0.7}, {95.5, 0.8}, {91.7, 1.0}, {84.4, 
            1.9}, {71.1, 3.9}, {58.9, 6.9}, {45.1, 10.8}, {37.1, 
            15.1}};

tie2=Interpolation[tbl4,InterpolationOrder -> 1]

tbl5 = 
      Table[Line[{{i, Extrc[i]}, {tie1[i], rafnt[tie1[i]]}}], {i, 1, 15, 
          2}]

Show[plt3, plt4, Graphics[{RGBColor[0, 1, 1], tbl5}], 
      Epilog -> Line[{{100, 0}, {0, 100}}], 
      PlotRange -> {{0, 100}, {0, 100}}]

Two Feeds and Solvent

F1 = 600. ; S = 1800. ; F2 = 400; xF1[1] = 65; xF1[2] = 35; xF1[3] = 0; 
    xF2[1] = 80; xF2[2] = 20; xF2[3] = 0;
xS[1] = 0; xS[2] = 0; xS[3] = 100. ;

Total Feed

sol0 = 
      Solve[{FT ==  F1 + F2, 
          FT xFT[1] ==  
            F1 xF1[1] + F2 xF2[1], 
FT xFT[2] ==  
            F1 xF1[2] + F2 xF2[2], 
          FT xFT[3] ==  F1 xF1[3] + F2 xF2[3]}, {FT,  xFT[1],  
          xFT[2],  xFT[3]}]

xFT[3] = sol0[[1, 1, 2]]

xFT[1] = sol0[[1, 2, 2]]

xFT[2] = sol0[[1, 3, 2]]

FT =  sol0[[1, 4, 2]]

Mixing Point

sol1 = 
      Solve[{M ==  FT + S, 
          M xM[1] ==  
            FT xFT[1] + S xS[1], 
M xM[2] ==  
            FT xFT[2] + S xS[2], 
          M xM[3] ==  FT xFT[3] + S xS[3]}, {M,  xM[1],  xM[2],  
          xM[3]}]

xM[1] = sol1[[1, 1, 2]]

xM[2] = sol1[[1, 2, 2]]

xM[3] = sol1[[1, 3, 2]]

M =  sol1[[1, 4, 2]]

Raffinate Product is such that acetic acid wt % is 10

xR[2] = 10

xR[1] = FindRoot[rafnt[x] ==  xR[2], {x, 80}][[1, 
          2]]

xR[3] = 100 - xR[1] - xR[2]

Extract Product

sys = {R + EP ==  M,  xR[1] R + xE[1] EP ==  M xM[1], 
        xR[2] R + xE[2] EP ==  M xM[2], 
        xR[3] R + xE[3] EP ==  M xM[3], 
        Extrc[xE[1]] ==  xE[2]}

sol2 = 
      FindRoot[sys, {EP, 100}, {R, 150}, {xE[1], 5}, {xE[2], 15}, {xE[3], 
          80}]

EP = sol2[[1, 2]]

R = sol2[[2, 2]]

xE[1] = sol2[[3, 2]]

xE[2] = sol2[[4, 2]]

xE[3] = sol2[[5, 2]]

Two Operating Points

sol3 = 
      Solve[{OP1 ==  FT - EP, 
          OP1 xOP1[1] ==  
            FT xFT[1] - 
              EP xE[1], 
OP1 xOP1[2] ==  
            FT xFT[2] - EP xE[2], 
          OP1 xOP1[3] ==  FT xFT[3] - EP xE[3]}, {OP1,  xOP1[1],  
          xOP1[2],  xOP1[3]}]

xOP1[1] = sol3[[1, 1, 2]]

xOP1[2] = sol3[[1, 2, 2]]

xOP1[3] = sol3[[1, 3, 2]]

OP1 = sol3[[1, 4, 2]]

sol3bis = 
      Solve[{OP2 ==  F2 - OP1, 
          OP2 xOP2[1] ==  
            F2 xF2[1] - 
              OP1 xOP1[1], 
OP2 xOP2[2] ==  
            F2 xF2[2] - OP1 xOP1[2], 
          OP2 xOP2[3] ==  F2 xF2[3] - OP1 xOP1[3]}, {OP2,  
          xOP2[1],  xOP2[2],  xOP2[3]}]

xOP2[1] = sol3bis[[1, 1, 2]]

xOP2[2] = sol3bis[[1, 2, 2]]

xOP2[3] = sol3bis[[1, 3, 2]]

OP2 = sol3bis[[1, 4, 2]]

Ploting Feed, Solvent, Extract, Raffinate and Operating Point

Show[plt3, plt4, 
      Graphics[{PointSize[0.02], RGBColor[1, 0, 0], Point[{xR[1], xR[2]}], 
          RGBColor[0, 1, 1], Point[{xM[1], xM[2]}], RGBColor[0, 0, 1], 
          Point[{xE[1], xE[2]}], RGBColor[0, 1, 1], 
          Point[{xOP1[1], xOP1[2]}], Point[{xOP2[1], xOP2[2]}], 
          RGBColor[0, 1, 1], Point[{xF1[1], xF1[2]}], 
          Point[{xF2[1], xF2[2]}], Point[{xFT[1], xFT[2]}], 
          Point[{xS[1], xS[2]}], 
          RGBColor[0, 1, 
            0], 
Line[{{xS[1], xS[2]}, {xFT[1], 
                xFT[2]}}], 
Line[{{xE[1], xE[2]}, {xR[1], 
                xR[2]}}], 
RGBColor[1, 0, 
            1], 
Line[{{xR[1], xR[2]}, {xOP1[1], 
                xOP1[2]}}], 
Line[{{xF2[1], 
                xF2[2]}, {xOP1[1], xOP1[2]}}], 
          Line[{{xFT[1], xFT[2]}, {xOP1[1], xOP1[2]}}]}], 
      Epilog -> Line[{{100, 0}, {0, 100}}], 
      PlotRange -> All]

Intersection between Line joining Operating Points and Extract Part of Equilibrium curve

OpLine[
        x_] := (x - 
              xOP1[1]) (xOP2[2] - xOP1[2])/(xOP2[1] - xOP1[1]) + 
        xOP1[2]

xinters = FindRoot[OpLine[x] ==  Extrc[x], {x, 1, 10}][[1, 
          2]]

Stepping off Stages

i = 1; xR[1, 1] = xR[1]; xR[1, 2] = xR[2]; xR[1, 3] = xR[3]; 
    xE[0, 2] = xE[2] - 10;

While[
      xE[i - 1, 2] <  xinters  &&  
        i < 100 , {
xE[i, 1] = 
          tie2[xR[i, 1]], 
xE[i, 2] = 
          Extrc[xE[i, 1]], 
xE[i, 3] = 
          100 - xE[i, 1] - 
            xE[i, 2], 
Sys = {El ==  Rl - OP1, 
            El xE[i, 1] ==  
              Rl xRl[1] - 
                OP1 xOP1[1], 
El xE[i, 2] == 
              Rl  xRl[2] - OP1 xOP1[2], 
            El xE[i, 3] ==  
              Rl xRl[3] - OP1 xOP1[3], 
xRl[2] ==  
              rafnt[xRl[1]]}, 
        sol[i] = 
          FindRoot[
            Sys, {El, 1700}, {Rl, 500}, {xRl[1], 80}, {xRl[2], 10}, {xRl[3], 
              5}], 
El[
            i] = sol[i][[1, 2]], 
Rl[
            i + 1] = sol[i][[2, 2]], 
xR[i + 1, 
            1] = sol[i][[3, 2]], 
xR[i + 1, 
            2] = sol[i][[4, 2]], 
xR[i + 1, 
            3] = sol[i][[5, 2]], i++}]

Nu1 = i;

While[
      xE[i - 1, 2] < xE[2]  &&  
        i < 100 , {
xE[i, 1] = 
          tie2[xR[i, 1]], 
xE[i, 2] = 
          Extrc[xE[i, 1]], 
xE[i, 3] = 
          100 - xE[i, 1] - 
            xE[i, 2], 
Sys = {El ==  Rl + OP2, 
            El xE[i, 1] ==  
              Rl xRl[1] + 
                OP2 xOP2[1], 
El xE[i, 2] == 
              Rl  xRl[2] + OP2 xOP2[2], 
            El xE[i, 3] ==  
              Rl xRl[3] + OP2 xOP2[3], 
xRl[2] ==  
              rafnt[xRl[1]]}, 
        sol[i] = 
          FindRoot[
            Sys, {El, 1700}, {Rl, 500}, {xRl[1], 60}, {xRl[2], 20}, {xRl[3], 
              5}], 
El[
            i] = sol[i][[1, 2]], 
Rl[
            i + 1] = sol[i][[2, 2]], 
xR[i + 1, 
            1] = sol[i][[3, 2]], 
xR[i + 1, 
            2] = sol[i][[4, 2]], 
xR[i + 1, 
            3] = sol[i][[5, 2]], i++}]

Ploting the Hunter and Nash Diagram. 2.89 equilibrium stages are needed to achieve the raffinate product requirement of 10 wt % acetic acid composition.

Nu2 = i - 1

Nu1

tbl5 = 
      Table[Point[{xE[i, 1], xE[i, 2]}], {i, 1, 
          Nu2}]
tbl6 = 
      Table[Point[{xR[i, 1], xR[i, 2]}], {i, 1, 
          Nu2}]
tbl7 = 
      Table[Line[{{xR[i, 1], xR[i, 2]}, {xOP1[1], xOP1[2]}}], {i, 1, 
          Nu1}]
tbl8 = 
      Table[Line[{{xE[i, 1], xE[i, 2]}, {xR[i, 1], xR[i, 2]}}], {i, 1, 
          Nu1}]

tbl9 = 
      Table[Line[{{xR[i, 1], xR[i, 2]}, {xOP2[1], xOP2[2]}}], {i, Nu1 + 1, 
          Nu2}]
tbl10 = 
      Table[Line[{{xE[i, 1], xE[i, 2]}, {xR[i, 1], xR[i, 2]}}], {i, Nu1 + 1, 
          Nu2}]

plt5 = 
      Show[plt3, plt4, 
        Graphics[{PointSize[0.02], RGBColor[0, 0, 1], Point[{xR[1], xR[2]}], 
            RGBColor[0, 1, 0], Point[{xE[1], xE[2]}], RGBColor[0, 0, 1], 
            Point[{xOP2[1], xOP2[2]}], Point[{xOP1[1], xOP1[2]}], 
            RGBColor[0, 0, 1], Point[{xF1[1], xF1[2]}], 
            Point[{xF2[1], xF2[2]}], Point[{xFT[1], xFT[2]}], 
            Point[{xS[1], xS[2]}], 
            RGBColor[0, 1, 0], 
tbl5, 
            tbl6, 
RGBColor[1, 0, 
              1], 
tbl7, 
            tbl9, 
RGBColor[0, 1, 
              0], 
tbl8, tbl10
}], 
        Epilog -> {Line[{{100, 0}, {0, 100}}], 
            Line[{{xOP1[1], xOP1[2]}, {xF2[1], xF2[2]}}]}, 
        AspectRatio -> 1.5, PlotRange -> {-10, 100}]

Show[plt5, PlotRange -> All]

(xE[2] - xE[2, 2])/(xE[3, 2] - xE[2, 2])

Cheking our results

S + F1 + F2

EP + R

Rl[2] + El[2]

El[1] + Rl[3] + F2