Applications of the Soave-Redlich-Kwong Equation of State
Authors' Data: Zakia Nasri and 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 Master's report of R. Privat under the direction of J-P Corriou.

θ =. ;

Off[NIntegrate::"nlim"]

Off[General::"spell1"]

Reduced Temperature

Tre = (θ + 273.15)/Tc;

Ethane's Critical Pressure, Critical Temperature and Accentric Factor. R is the Universal Gas Constant

Tc = 273.15 + 32.28;  Pc = 48.84;  ω = 0.09860; 
    R = 83.14;

Soave-Redlich-Kwong Equation of State

P[T_, v_] := R T/(v - b) - a/(v (v + b))

b = 0.08664 R Tc/Pc

m = 0.480 + 1.574 ω - 0.176 ω^2

a = 
      0.42748 (R Tc)^2/Pc (1 + m (1 - (Sqrt[Tre])))^2

tbl = 
        Table[P[θ + 273.15, v], {θ, 1, 30, 5}];

Plotting Isotherms

plt[θ_] := 
      Plot[tbl[[θ]], {v, 0, 2500}, 
        PlotRange -> {{0, 2500}, {0, 100}}, 
        PlotStyle -> Hue[θ/6], 
        DisplayFunction -> Identity]

plot1 = 
      Show[Table[plt[θ], {θ, 1, 6}], 
        DisplayFunction -> $DisplayFunction]

Finding vapor pressure and liquid and gas molar volumes at different temperatures. Dots in the figure show vapor liquid equilibrium curve for ethane.

sol[θ_] := 
        FindRoot[{NIntegrate[-P[273.15 + θ, v] + 
                  P[273.15 + θ, x], {v, x, z}] == 
              NIntegrate[
                P[273.15 + θ, v] - P[273.15 + θ, y], {v, z, 
                  y}] , P[273.15 + θ, x] == 
              P[273.15 + θ, z], 
            P[273.15 + θ, z] == P[273.15 + θ, y]}, {x, 
            50}, {z, 120}, {y, 400}, MaxIterations -> 100]

tbl1 = 
      Table[{x /. sol[θ], 
          P[θ + 273.15, x /. sol[θ]]}, {θ, 1, 30, 
          5}]

tbl2 = 
      Table[{y /. sol[θ], 
          P[θ + 273.15, y /. sol[θ]]}, {θ, 1, 30, 
          5}]

TBL = Join[tbl1, Sort[tbl2]]

plot2 = ListPlot[TBL, PlotStyle -> PointSize[0.025]]

Show[plot1, plot2, PlotRange -> {{0, 600}, {5, 50}}]

The equilibrium data is obtained by solving a system of equation so that the two areas are equal.

θ = 11;

plt5 = 
      Show[plt[3], 
        Epilog -> {Line[{{x /. sol[11], 
                  P[11 + 273.15, x /. sol[11]]}, {y /. sol[11], 
                  P[11 + 273.15, y /. sol[11]]}}], {PointSize[0.025], 
              Point[{y /. sol[11], P[11 + 273.15, x /. sol[11]]}], 
              Point[{x /. sol[11], P[11 + 273.15, x /. sol[11]]}], 
              Point[{z /. sol[11], P[11 + 273.15, x /. sol[11]]}]}}, 
        DisplayFunction -> $DisplayFunction, 
        PlotRange -> {{0, 600}, {5, 50}}]

NIntegrate[(-P[273.15 + θ, v] + 
              P[273.15 + θ, x]) /. θ -> 11, {v, x, z}] /. 
      sol[11]

NIntegrate[(P[273.15 + θ, v] - 
              P[273.15 + θ, y]) /. θ -> 11, {v, z, y}] /. 
      sol[11]

<< Graphics`FilledPlot`

plt7 = 
      FilledPlot[{P[θ + 273.15, v] /. θ -> 
              11, P[273.15 + θ, y] /. θ -> 11 /. 
            sol[11]}, {v, x, z} /. sol[11], 
        DisplayFunction -> Identity]

plt8 = 
      FilledPlot[{P[θ + 273.15, v] /. θ -> 
              11, P[273.15 + θ, y] /. θ -> 11 /. 
            sol[11]}, {v, z, y} /. sol[11], 
        DisplayFunction -> Identity]

Show[plt5, plt7, plt8]

Vapor pressure data versus temperature. Comparison with experimentally derived expressions of vapor pressure. We get a perfect fit with the modified Antoine's equation.

plot3 = 
      ListPlot[Table[{θ, 
            P[θ + 273.15, x /. sol[θ]]}, {θ, 1, 30, 5}], 
        PlotStyle -> PointSize[0.025]]

Modified Antoine's equation obtained from HYSYS 3.2

plot4 = 
      Plot[Exp[(44.0103 - 2.56882 10^3/(θ + 273.15) - 
                4.97635 Log[(θ + 273.15)] + 
                1.46447 10^-5 (θ + 273.15)^2)]/
          100, {θ, 0, 30}, 
        PlotStyle -> RGBColor[0, 1, 0]]

Antoine's equation obtained from Equilibrium-Stage Separation Operations in Chemical Engineering by Henley and Seader, Wiley, 1981.

plt5 = 
      Plot[Pc Exp[
            5.383894 - 
              2847.921/(434.898 + 1.8 θ + 32)], {θ, 0, 
          30}, PlotStyle -> RGBColor[0, 0, 1]]

Show[plot3, plt5, plot4, PlotRange -> {20, 60}]