Liquid Liquid Equilibrium Using NDSolve
Case of a ternary system: water - ethyl acetate - acetone at 15 °C
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Off[General::"spell1"]

water = component 1
ethylacetate = component 2
acetone = component 3

Activities using the NRTL model

GAMMA1 = 
      Exp[(t21/(R T) G21 x2 + t31/(R T) G31 x3)/(x1 + 
                G21 x2 + G31 x3) + 
          x1/(x1 + G21 x2 + 
                  G31 x3) (-(x2 t21/(R T) G21 + 
                      x3 t31/(R T) G31)/(x1 + G21 x2 + 
                    G31 x3)) + 
          x2 G12/(G12 x1 + x2 + 
                  G32 x3) (t12/(R T) - (x1 t12/(R T) G12 + x3 t32/(R T) G32)/(G12 x1 + x2 + 
                      G32 x3)) + 
x3 G13/(G13 x1 + 
                  G23 x2  + 
                  x3) (t13/(R T) - (x1 t13/(R T) G13 + 
                      x2 t23/(R T) G23)/(G13 x1 + G23 x2  + 
                      x3))]

GAMMA1p = 
      Exp[(t21/(R T) G21 xp2 + t31/(R T) G31 xp3)/(xp1 + 
                G21 xp2 + G31 xp3) + 
          xp1/(xp1 + G21 xp2 + 
                  G31 xp3) (-(xp2 t21/(R T) G21 + 
                      xp3 t31/(R T) G31)/(xp1 + G21 xp2 + 
                    G31 xp3)) + 
          xp2 G12/(G12 xp1 + xp2 + 
                  G32 xp3) (t12/(R T) - (xp1 t12/(R T) G12 + xp3 t32/(R T) G32)/(G12 xp1 + xp2 + 
                      G32 xp3)) + 
xp3 G13/(G13 xp1 + G23 xp2  + 
                  xp3) (t13/(R T) - (xp1 t13/(R T) G13 + 
                      xp2 t23/(R T) G23)/(G13 xp1 + G23 xp2  + 
                      xp3))]

GAMMA2 = 
      Exp[(t12/(R T) G12 x1 + 
                t32/(R T) G32 x3)/(G12 x1 + x2 + G32 x3) + 
          x1 G21/(x1 + G21 x2 + 
                  G31 x3) (t21/(R T) - (x2 t21/(R T) G21 + x3 t31/(R T) G31)/(x1 + G21 x2 + G31 x3)) + 
          x2/(G12 x1 + x2 + 
                  G32 x3) (-(x1 t12/(R T) G12 + 
                      x3 t32/(R T) G32)/(G12 x1 + x2 + 
                    G32 x3)) + 
          x3 G23/(G13 x1 + G23 x2  + 
                  x3) (t23/(R T) - (x1 t13/(R T) G13 + 
                      x2 t23/(R T) G23)/(G13 x1 + G23 x2  + 
                      x3))]

GAMMA2p = 
      Exp[(t12/(R T) G12 xp1 + 
                t32/(R T) G32 xp3)/(G12 xp1 + xp2 + G32 xp3) + 
          xp1 G21/(xp1 + G21 xp2 + 
                  G31 xp3) (t21/(R T) - (xp2 t21/(R T) G21 + xp3 t31/(R T) G31)/(xp1 + G21 xp2 + G31 xp3)) + 
          xp2/(G12 xp1 + xp2 + 
                  G32 xp3) (-(xp1 t12/(R T) G12 + 
                      xp3 t32/(R T) G32)/(G12 xp1 + xp2 + 
                    G32 xp3)) + 
          xp3 G23/(G13 xp1 + G23 xp2  + 
                  xp3) (t23/(R T) - (xp1 t13/(R T) G13 + 
                      xp2 t23/(R T) G23)/(G13 xp1 + G23 xp2  + 
                      xp3))]

GAMMA3 = 
      Exp[(t13/(R T) G13 x1 + 
                t23/(R T) G23 x2)/(G13 x1 + G23 x2 + x3) + 
          x1 G31/(x1 + G21 x2 + 
                  G31 x3) (t31/(R T) - (x2 t21/(R T) G21 + x3 t31/(R T) G31)/(x1 + G21 x2 + G31 x3)) + 
          x2 G32/(G12 x1 + x2 + 
                  G32 x3) (t32/(R T) - (x1 t12/(R T) G12 + x3 t32/(R T) G32)/(G12 x1 + x2 + G32 x3)) + 
          x3/(G13 x1 + G23 x2  + 
                  x3) (-(x1 t13/(R T) G13 + 
                      x2 t23/(R T) G23)/(G13 x1 + G23 x2  + 
                    x3))]

GAMMA3p = 
      Exp[(t13/(R T) G13 xp1 + 
                t23/(R T) G23 xp2)/(G13 xp1 + G23 xp2 + xp3) + 
          xp1 G31/(xp1 + G21 xp2 + 
                  G31 xp3) (t31/(R T) - (xp2 t21/(R T) G21 + xp3 t31/(R T) G31)/(xp1 + G21 xp2 + G31 xp3)) + 
          xp2 G32/(G12 xp1 + xp2 + 
                  G32 xp3) (t32/(R T) - (xp1 t12/(R T) G12 + xp3 t32/(R T) G32)/(G12 xp1 + xp2 + G32 xp3)) + 
          xp3/(G13 xp1 + G23 xp2  + 
                  xp3) (-(xp1 t13/(R T) G13 + 
                      xp2 t23/(R T) G23)/(G13 xp1 + 
                    G23 xp2  + xp3))]

G12 = Exp[-a12 t12/(R T)]; 
    G21 = Exp[-a12 t21/(R T)]; G13 = Exp[-a13 t13/(R T)]; 
    G31 = Exp[-a13 t31/(R T)]; G32 = Exp[-a32 t32/(R T)]; 
    G23 = Exp[-a32 t23/(R T)];

Gas constant, temperature and binary interaction parameters for the NRTL model

R = 1.987;
a13 = 0.585600018501282; 
    a12 = 0.410400003194809;
a32 = 0.291200011968613;
t31 = 750.318115234375; 
    t13 = 1299.39501953125;
t23 = 467.125701904297; 
    t32 = -304.988403320313;
t12 = 2316.36303710938;
t21 = 935.687927246094;

T = 273.15 + 15

Computing a Tie line with FindRoot

X1 = 0.7; X2 = 0.2;

\(sol[1] = 
      FindRoot[{\((GAMMA1\ x1 == \ 
              GAMMA1p\ xp1)\), \((GAMMA2\ x2 == \ 
              GAMMA2p\ xp2)\), 
\((GAMMA3\ x3 == \ 
              GAMMA3p\ xp3)\), 
\((x1\ L + 
                xp1\ Lp == \ X1)\), \((x2\ L + xp2\ Lp == \ 
              X2)\), \((1 == \ xp1 + xp2 + xp3)\), \((1 == \ 
              x1 + x2 + x3)\), \((L + Lp == 1)\)}, {x1, 0.95}, {x2, 
          0.01}, {xp1, 0.5}, {xp2, 0.4}, {L, 0.05}, {Lp, 0.04}, {x3, 
          0.01}, {xp3, 0.1}, 
	MaxIterations -> 10000]\)

Computing a Tie line with NDSolve

Coordinates of the first point on the tie line.

x1 = sol[1][[1, 2]]

x2 = sol[1][[2, 2]]

x3 = 1 - x1 - x2

Using NDSolve to look for the steady state, which corresponds to the solution of the nonlinear algebraic equations.

sol1 = 
      NDSolve[{D[
              xp1[t], {t, 
                1}] ==  (GAMMA1 x1 -  GAMMA1p xp1 /. 
                    xp1 -> xp1[t] /. xp2 ->  xp2[t] /. 
                xp3 ->  xp3[t]), 
D[
              xp2[t], {t, 
                1}] ==  (GAMMA2 x2 -  GAMMA2p xp2 /. 
                    xp1 -> xp1[t] /. xp2 ->  xp2[t] /. 
                xp3 ->  xp3[t]), 
D[
              xp3[t], {t, 
                1}] ==  (GAMMA3 x3 -  GAMMA3p xp3 /. 
                    xp1 -> xp1[t] /. xp2 ->  xp2[t] /. 
                xp3 ->  xp3[t]), 
xp1[0] ==  
            0.4, xp2[0] ==  0.2, xp3[0] ==  0.1}, {xp1, xp2, 
          xp3}, {t, 0, 1000}]

Coordinates of the second point on the tie line are the steady state solution (obtained at t=1000). We get the same result as the one found using FindRoot (xp1, xp2, xp3)

First[Evaluate[xp1[t] /. sol1 /. t ->  1000]]

First[Evaluate[xp2[t] /. sol1 /. t ->  1000]]

First[Evaluate[xp3[t] /. sol1 /. t ->  1000]]

We reach the steady state rapidly (we are almost at the steady state at t=4).

Plot[{Evaluate[xp1[t] /. sol1], Evaluate[xp2[t] /. sol1], 
        Evaluate[xp3[t] /. sol1]}, {t, 0, 4}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0], 
          RGBColor[0, 0, 1]}]