Vapor Phase Solubility of Decane in Nitrogen at 50°C
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :
This problem was presented in the excellent book by Prausnitz et al., Molecular Thermodynamics of Fluid-Phase Equilibria, 3rd edition, Prentice Hall, 1999.

Off[General::"spell1"]

Reduced Pressure and Temperature

Tre[i_] := T/Tc[i]

Pre[i_] := P/Pc[i]

Equilibrium constants using the ϕ-ϕ methode and the Peng Robinson Equation of State. Binary interaction parameters for the Peng Robinson EOS were taken from HYSYS 3.2.

ap[i_] := 0.45724 a[i]  R^2 Tc[i]^2/Pc[i]

Ap[i_] := ap[i] P/(R T)^2;

bp[i_] := 0.07780 R Tc[i]/Pc[i]

Bp[i_] := P bp[i]/(R T)

(Subscript[ϕ, v])[i_] := 
      Exp[((Subscript[Z, v]) - 1) Bp[i]/(Subscript[B, v]) - Log[(Subscript[Z, v]) - (Subscript[B, v])] - 
          (Subscript[A, v])/(2 (Sqrt[2]) 
                  (Subscript[B, v])) (2 (∑_(j = 1)^2) y[j] Ab[i, j]/
                        (Subscript[A, v]) - 
                Bp[i]/(Subscript[B, v])) Log[((Subscript[Z, v]) + (1 + (Sqrt[2])) 
                      (Subscript[B, v]))/((Subscript[Z, v]) + (1 - (Sqrt[2])) (Subscript[B, v]))]]

(Subscript[ϕ, l])[i_] := 
      Exp[((Subscript[Z, l]) - 1) Bp[i]/(Subscript[B, l]) - Log[(Subscript[Z, l]) - (Subscript[B, l])] - 
          (Subscript[A, l])/(2 (Sqrt[2]) 
                  (Subscript[B, l])) (2 (∑_(j = 1)^2) x[j] Ab[i, j]/
                        (Subscript[A, l]) - 
                Bp[i]/(Subscript[B, l])) Log[((Subscript[Z, l]) + (1 + (Sqrt[2])) 
                      (Subscript[B, l]))/((Subscript[Z, l]) + (1 - (Sqrt[2])) (Subscript[B, l]))]]

a[i_] := (1 + m[i] (1 - Tre[i]^0.5))^2

m[i_] := 
      0.37464 + 1.54226 ω[i] - 0.26992 ω[i]^2

(Subscript[A, v]) = (∑_(i = 1)^2)(∑_(j = 1)^2) y[i] y[j] Ab[i, 
                j];

Ab[i_, j_] = (Ap[i] Ap[j])^0.5 (1 - k[i, j]);

(Subscript[B, v]) = (∑_(i = 1)^2) y[i] Bp[i];

(Subscript[A, l]) = (∑_(i = 1)^2)(∑_(j = 1)^2) x[i] x[j] Ab[i, 
                j];

(Subscript[B, l]) = (∑_(i = 1)^2) x[i] Bp[i];

K[i_] := (Subscript[ϕ, l])[i]/(Subscript[ϕ, v])[i];

k[1, 1] = 0; k[2, 2] = 0;
k[2, 1] = k[1, 2];
k[1, 2] = 0.1316;

R = 1.987;

Critical Pressure and Temperature and acentric factor for n-Decane and Nitrogen

Pc[1] = 306.0; Tc[1] = 1114.0; ω[1] = 0.4869;

Pc[2] = 492.9; Tc[2] = 227.3; ω[2] = 0.0206;

Flash Pressure and Temperature are  77°F and P is varied from 0.5 to 100 bars. 
Feed is arbitrarily composed of 50% n-Decane and 50% Nitrogen.

T = 25 1.8 + 32 + 459.67

z[1] := 0.5; z[2] = 0.5;

Rachford and Rice Equation and Flash Calculations.

Eq[i_] := y[i] ==  K[i] x[i]

V = (Subscript[ϕ, vap]) F; L = (1 - (Subscript[ϕ, vap])) F;

EQ[i_] := 
        x[i] ==  
          z[i]/(1 + (Subscript[ϕ, vap]) (K[i] - 1));

sol[0] = {x[2] -> 0.0001, 
        x[1] -> 0.9999, (Subscript[ϕ, vap]) -> 0.8, y[2] -> 0.999, 
        y[1] -> 0.001, (Subscript[Z, l]) -> 0.01, (Subscript[Z, v]) -> 0.99}

For[i = 1, 
      i < 17, {P = 
          14.696 {0.5, 1, 2, 4, 6, 8, 10, 20, 30, 40, 50, 60, 70, 80, 90, 
                100}[[i]], 
        sol[i] = 
          FindRoot[{Eq[1], (∑_(i = 1)^2) y[i] == 1, 
                EQ[1], (∑_(i = 1)^2) x[i] == 
                  1, (∑_(i = 1)^2)(( z[i] (K[i] - 1))/(1 + (Subscript[ϕ, vap]) (K[i] - 1))) == 0, 
                (Subscript[Z, l])^3 + (-1 + (Subscript[B, l])) (Subscript[Z, l])^2 + 
                    (Subscript[Z, l]) ((Subscript[A, l]) - 3  (Subscript[B, l])^2 - 2  (Subscript[B, l])) - (Subscript[A, l]) (Subscript[B, l]) + 
                    (Subscript[B, l])^2 + (Subscript[B, l])^3 == 0, 
                (Subscript[Z, v])^3 + (-1 + (Subscript[B, v])) (Subscript[Z, v])^2 + 
                    (Subscript[Z, v]) ((Subscript[A, v]) - 3  (Subscript[B, v])^2 - 2  (Subscript[B, v])) - (Subscript[A, v]) (Subscript[B, v]) + 
                    (Subscript[B, v])^2 + (Subscript[B, v])^3 == 0}, {x[1], 
                x[1] /. sol[i - 1]}, {x[2], 
                x[2] /. sol[i - 1]}, {(Subscript[ϕ, vap]), (Subscript[ϕ, vap]) /. 
                  sol[i - 1]}, {y[1], y[1] /. sol[i - 1]}, {y[2], 
                y[2] /. sol[i - 1]}, {(Subscript[Z, l]), (Subscript[Z, l]) /. sol[i - 1]}, {(Subscript[Z, v]), 
                (Subscript[Z, v]) /. sol[i - 1]}, MaxIterations -> 10000] // Chop, 
        Print[i], i++}]

P = 
      14.696 {0.5, 1, 2, 4, 6, 8, 10, 20, 30, 40, 50, 60, 70, 80, 90, 
          100}

Y = Table[y[1] /. sol[i], {i, 1, 16}]

tbl = Transpose[{P, Y}]

<< Graphics`Graphics`

LogListPlot[tbl, PlotStyle -> RGBColor[1, 0, 0], 
      PlotJoined -> True, PlotRange -> {0.0001, 0.01}, 
      AspectRatio -> 1.15]

Conclusion: Results show only qualitative agreement with Prausnitz et al. (1999). The binary interaction parameter, taken from HYSYS 3.2, may be responsible of this small discrepancy.