Calculation of Solubity of Naphtalene in Supercritical CO2 using 
the Peng Robinson EOS
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 S. I. Sandler, Chemical and Engineering Thermodynamics, 3rd Edition, Wiley 1999 (Illustration 8.5-6 page 583).

Off[General::"spell1"]

Reduced Pressure and Temperature

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

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

Fugacity coefficient are obtained using 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]))]]

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

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

R = 1.987;

Critical Pressure and Temperature and acentric factor for Naphtalene and CO2

Pc[1] = 576.1; Tc[1] = 1347.0; ω[1] = 0.2934;

Pc[2] = 1070.5; Tc[2] = 547.6; ω[2] = 0.177;

Vapor pressure for Naphtalene and CO2

Psat[2] = 
      Pc[2] Exp[Asat[2] - Bsat[2]/(T - 459.67 + Csat[2])]

Asat[2] = 6.470232; Bsat[2] = 3521.259; Csat[2] = 455.869;

Psat[1] = 14.696 10^(8.722 - 3783/(T/1.8))

Solid molar volume

VS = 125 1/28317.  453.593

Temperature is 35°C

T = 32 + 35 1.8 + 459.67

Solubility of Naphtalene in CO2 at various pressures 
( P=1275.61, 1443.15, 1954.57, 2931.85 and 3751.89 psia )

P =  86.8 14.696

S1 = 
      FindMinimum[(y[1] (P (Subscript[ϕ, v])[1]) -  
                Psat[1] Exp[
                    VS (P - 
                            Psat[1])/(10.73 T)])^2 + ((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)^2 + (1 - y[1] - 
                y[2])^2, {y[1], 0.006}, {y[2], 0.9}, {(Subscript[Z, v]), 0.85}, 
        MaxIterations ->  100000]

P =  98.2 14.696

S2 = 
      FindMinimum[(y[1] (P (Subscript[ϕ, v])[1]) -  
                Psat[1] Exp[
                    VS (P - 
                            Psat[1])/(10.73 T)])^2 + ((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)^2 + (1 - y[1] - 
                y[2])^2, {y[1], 0.006}, {y[2], 0.9}, {(Subscript[Z, v]), 0.85}, 
        MaxIterations ->  100000]

P =  133.0 14.696

S3 = 
      FindMinimum[(y[1] (P (Subscript[ϕ, v])[1]) -  
                Psat[1] Exp[
                    VS (P - 
                            Psat[1])/(10.73 T)])^2 + ((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)^2 + (1 - y[1] - 
                y[2])^2, {y[1], 0.006}, {y[2], 0.9}, {(Subscript[Z, v]), 0.85}, 
        MaxIterations ->  100000]

P =  199.5 14.696

S4 = 
      FindMinimum[(y[1] (P (Subscript[ϕ, v])[1]) -  
                Psat[1] Exp[
                    VS (P - 
                            Psat[1])/(10.73 T)])^2 + ((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)^2 + (1 - y[1] - 
                y[2])^2, {y[1], 0.02}, {y[2], 0.9}, {(Subscript[Z, v]), 0.45}, 
        MaxIterations ->  100000]

P =  255.3 14.696

S5 = 
      FindMinimum[(y[1] (P (Subscript[ϕ, v])[1]) -  
                Psat[1] Exp[
                    VS (P - 
                            Psat[1])/(10.73 T)])^2 + ((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)^2 + (1 - y[1] - 
                y[2])^2, {y[1], 0.02}, {y[2], 0.9}, {(Subscript[Z, v]), 0.85}, 
        MaxIterations ->  100000]

Ideal solubility of Naphtalene in CO2 at various pressures 
( P=1275.61, 1443.15, 1954.57, 2931.85 and 3751.89 psia )

P = 14.696 {86.8, 98.2, 133.0, 199.5, 255.3}

Psat[1]/P

Enhancement factor at various pressures

{y[1] /. S1[[2]], y[1] /. S2[[2]], y[1] /. S3[[2]], 
        y[1] /. S4[[2]], y[1] /. S5[[2]]}/%

Poynting factor at various pressures

Exp[VS (P - Psat[1])/(10.73 T)]

Plot of predicted and experimental solubility of Naphtalene in CO2 versus various pressure.

<< Graphics`Graphics`

tbl = 
      Transpose[{14.696 {86.8, 98.2, 133.0, 199.5, 255.3}, {y[1] /. 
              S1[[2]], y[1] /. S2[[2]], y[1] /. S3[[2]], 
            y[1] /. S4[[2]], y[1] /. S5[[2]]}}]

plt1 = 
      LogListPlot[tbl, PlotJoined ->  True, 
        PlotStyle ->  RGBColor[0, 0, 1], 
        DisplayFunction -> Identity]

tblexp = 
      Transpose[{14.696 {86.8, 98.2, 133.0, 199.5, 255.3}, {0.00750, 
            0.00975, 0.01410, 0.01709, 0.01922}}]

plt2 = 
      LogListPlot[tblexp, 
        PlotStyle ->  {PointSize[0.02], RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]

Show[plt1, plt2, PlotRange -> All, 
      DisplayFunction ->  $DisplayFunction]

Temperature is 60.4°C

T = 32 + 60.4 1.8 + 459.67

Solubility of Naphtalene in CO2 at various pressures 
( P=1593.05, 1966.32, 2241.14, 2413.08, 2830.45 and 3027.38 psia )

P =  108.4 14.696

S1 = 
      FindMinimum[(y[1] (P (Subscript[ϕ, v])[1]) -  
                Psat[1] Exp[
                    VS (P - 
                            Psat[1])/(10.73 T)])^2 + ((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)^2 + (1 - y[1] - 
                y[2])^2, {y[1], 0.006}, {y[2], 0.9}, {(Subscript[Z, v]), 0.85}, 
        MaxIterations ->  100000]

P =  133.8 14.696

S2 = 
      FindMinimum[(y[1] (P (Subscript[ϕ, v])[1]) -  
                Psat[1] Exp[
                    VS (P - 
                            Psat[1])/(10.73 T)])^2 + ((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)^2 + (1 - y[1] - 
                y[2])^2, {y[1], 0.006}, {y[2], 0.9}, {(Subscript[Z, v]), 0.85}, 
        MaxIterations ->  100000]

P =  152.5 14.696

S3 = 
      FindMinimum[(y[1] (P (Subscript[ϕ, v])[1]) -  
                Psat[1] Exp[
                    VS (P - 
                            Psat[1])/(10.73 T)])^2 + ((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)^2 + (1 - y[1] - 
                y[2])^2, {y[1], 0.006}, {y[2], 0.9}, {(Subscript[Z, v]), 0.85}, 
        MaxIterations ->  100000]

P =  164.2 14.696

S4 = 
      FindMinimum[(y[1] (P (Subscript[ϕ, v])[1]) -  
                Psat[1] Exp[
                    VS (P - 
                            Psat[1])/(10.73 T)])^2 + ((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)^2 + (1 - y[1] - 
                y[2])^2, {y[1], 0.04}, {y[2], 0.9}, {(Subscript[Z, v]), 0.45}, 
        MaxIterations ->  100000]

P =  192.6 14.696

S5 = 
      FindMinimum[(y[1] (P (Subscript[ϕ, v])[1]) -  
                Psat[1] Exp[
                    VS (P - 
                            Psat[1])/(10.73 T)])^2 + ((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)^2 + (1 - y[1] - 
                y[2])^2, {y[1], 0.055}, {y[2], 0.9}, {(Subscript[Z, v]), 0.85}, 
        MaxIterations ->  100000]

P =  206.0 14.696

S6 = 
      FindMinimum[(y[1] (P (Subscript[ϕ, v])[1]) -  
                Psat[1] Exp[
                    VS (P - 
                            Psat[1])/(10.73 T)])^2 + ((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)^2 + (1 - y[1] - 
                y[2])^2, {y[1], 0.02}, {y[2], 0.9}, {(Subscript[Z, v]), 0.85}, 
        MaxIterations ->  100000]

Ideal solubility of Naphtalene in CO2 at various pressures 
( P=1593.05, 1966.32, 2241.14, 2413.08, 2830.45 and 3027.38 psia )

P = 14.696 {108.4, 133.8, 152.5, 164.2, 192.6, 206.0}

Psat[1]/P

Enhancement factor at various pressures

{y[1] /. S1[[2]], y[1] /. S2[[2]], y[1] /. S3[[2]], 
        y[1] /. S4[[2]], y[1] /. S5[[2]], 
        y[1] /. S6[[2]]}/%

Poynting factor at various pressures

Exp[VS (P - Psat[1])/(10.73 T)]

Plot of predicted and experimental solubility of Naphtalene in CO2 versus various pressure.

<< Graphics`Graphics`

tbl = 
      Transpose[{14.696 {108.4, 133.8, 152.5, 164.2, 192.6, 206.0}, {y[1] /. 
              S1[[2]], y[1] /. S2[[2]], y[1] /. S3[[2]], 
            y[1] /. S4[[2]], y[1] /. S5[[2]], 
            y[1] /. S6[[2]]}}]

plt1 = 
      LogListPlot[tbl, PlotJoined ->  True, 
        PlotStyle ->  RGBColor[0, 0, 1], 
        DisplayFunction -> Identity]

tblexp = 
      Transpose[{14.696 {108.4, 133.8, 152.5, 164.2, 192.6, 
              206.0}, {0.00524, 0.01516, 0.02589, 0.04296, 0.05386, 
            0.06259}}]

plt2 = 
      LogListPlot[tblexp, 
        PlotStyle ->  {PointSize[0.02], RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity]

Show[plt1, plt2, PlotRange -> All, 
      DisplayFunction ->  $DisplayFunction]