Study of a Binary Mixture of Methane and Propane
using the Virial Equation of State
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 Thermodynamique, Application au Génie Chimique et à l'Industrie Pétrolière by Jean Vidal, Technip, 1997.

Off[InterpolatingFunction::"dmval"]

Off[Power::"infy"]

Off[∞::"indet"]

Experimental molar volume in (cm^3)/mol as a function of composition for binary mixture Methane-Propane at T=344.15 K and P=1.377 MPa. R is the universal gas constant.

P = 1.377;

T = 344.15;

R = 8.314;

vexp = {1748, 1791, 1835, 1873, 1910, 1941, 1967, 1992, 2016, 2035, 
          2047};

y = Table[i, {i, 0, 1, 0.1}]

Join[{{"y", "vexp"}}, Transpose[{y, vexp}]] // 
      TableForm

Determination of second virial coefficients ((Subscript[B, 11]), (Subscript[B, 12]) and (Subscript[B, 22])) using experimental molar volumes and the least squares method.

vtheo = 
      Map[R T/P + B11 #^2 + 2 B12 # (1.  - #) + 
            B22 (1 - #)^2 &, y]

sol1 = 
      FindMinimum[
        Sum[(vexp[[i]] - vtheo[[i]])^2, {i, 1, 
            11}], {B11, -5}, {B12, -200}, {B22, -240}]

Plotting theoretically determined and experimental molar volumes.

plt1 = 
      ListPlot[Transpose[{y, vexp}], 
        PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.02]}, 
        DisplayFunction -> Identity]

plt2 = 
      ListPlot[Transpose[{y, vtheo /. sol1[[2]]}], 
        PlotStyle -> RGBColor[0, 0, 1], PlotJoined -> True, 
        DisplayFunction -> Identity]

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

Computing partial molar volumes of Methane and Propane as a function of composition (expressed in (cm^3)/mol).

vpart1 = 
      Map[R T/P + 
              2 (B11 # + B12 (1 - #)) - (B11 #^2 + 
                  2 B12 # (1.  - #) + B22 (1 - #)^2) &, y] /. 
        sol1[[2]]

vpart2 = 
      Map[R T/P + 
              2 (B12 # + B22 (1 - #)) - (B11 #^2 + 
                  2 B12 # (1.  - #) + B22 (1 - #)^2) &, y] /. 
        sol1[[2]]

Plotting partial molar volumes.

plt3 = 
      ListPlot[Transpose[{y, vpart1}], PlotJoined -> True, 
        PlotStyle -> {RGBColor[0, 1, 0], PointSize[0.02]}, 
        DisplayFunction -> Identity]

plt4 = 
      ListPlot[Transpose[{y, vpart2}], PlotJoined -> True, 
        PlotStyle -> {RGBColor[0, 1, 1], PointSize[0.02]}, 
        DisplayFunction -> Identity]

Show[plt3, plt4, DisplayFunction -> $DisplayFunction]

Show[plt1, plt2, plt3, plt4, 
      DisplayFunction -> $DisplayFunction]

Computing fugacity coefficients and fugacities (expressed in MPa) of Methane and Propane as a function of composition.

(Subscript[φ, 1]) = 
      Map[Exp[1/(R T) (2 (B11 # + 
                          B12 (1 - #)) - (B11 #^2 + 
                        2 B12 # (1.  - #) + 
                        B22 (1 - #)^2))] &, y] /. 
        sol1[[2]]

(Subscript[f, 1]) = P y (Subscript[φ, 1])

(Subscript[φ, 2]) = 
      Map[Exp[1/(R T) (2 (B12 # + 
                          B22 (1 - #)) - (B11 #^2 + 
                        2 B12 # (1.  - #) + 
                        B22 (1 - #)^2))] &, y] /. 
        sol1[[2]]

(Subscript[f, 2]) = P (1 - y) (Subscript[φ, 2])

\(Join[{{"y", \*""(Subscript[φ, 1])"", \*""(Subscript[φ, 2])"", \*""(Subscript[f, 1])"", \*""(Subscript[f, 2])""}}, Transpose[{y, φ\_1, φ\_2, f\_1, f\_2}]] // 
      TableForm\)

Fugacities of pure components at reference state in MPa

\(\(f\_1\^*\) = f\_1[\([11]\)]\)

\(\(f\_2\^*\) = f\_2[\([1]\)]\)

Activities and activity coefficients

\(a\_1 = f\_1/\(f\_1\^*\)\)

\(a\_2 = f\_2/\(f\_2\^*\)\)

(Subscript[γ, 1]) = (Subscript[a, 1])/y

(Subscript[γ, 2]) = (Subscript[a, 2])/(1 - y)

gamma1 = 
      Interpolation[
        Table[{y[[i]], (Subscript[γ, 1])[[i]]}, {i, 2, 
            11}], 
InterpolationOrder  -> 1]

gamma1[0]

(Subscript[γ, 1]) = Join[{gamma1[0]}, Drop[(Subscript[γ, 1]), 1]]

gamma2 = 
      Interpolation[
        Table[{1 - y[[i]], (Subscript[γ, 2])[[i]]}, {i, 1, 10}], 
        InterpolationOrder  -> 1]

gamma2[0]

(Subscript[γ, 2]) = Join[Drop[(Subscript[γ, 2]), -1], {gamma2[0]}]

\(Join[{{"y", \*""(Subscript[a, 1])"", \*""(Subscript[a, 2])"", \*""(Subscript[γ, 1])"", \*""(Subscript[γ, 2])""}}, 
        Transpose[{y, a\_1, a\_2, γ\_1, γ\_2}]] // 
      TableForm\)

Computing and plotting excess gibbs free energy. Pure component fugacity coefficients, for Methane and Propane, are \(\(φ\_1\^*\)\) and \(\(φ\_2\^*\)\) .

\(\(φ\_1\^*\) = \(f\_1\^*\)/P\)

\(\(φ\_2\^*\) = \(f\_2\^*\)/P\)

\(gex = 
      R\ T \((y\ Log[φ\_1] + \((1 - 
                    y)\)\ Log[φ\_2] - \((y\ Log[\(φ\_1\^*\)] + \((1 - y)\)\ Log[\(φ\_2\^*\)])\))\) // Chop\)

plt5 = 
      ListPlot[Transpose[{y, gex}], PlotStyle -> RGBColor[1, 0, 0], 
        PlotJoined -> True]

Computing and plotting excess chemical potentials

(Subscript[μ, 1  ex]) = R T Log[(Subscript[γ, 1])]

(Subscript[μ, 2  ex]) = R T Log[(Subscript[γ, 2])]

plt6 = 
      ListPlot[Transpose[{y, (Subscript[μ, 1  ex])}], PlotJoined -> True, 
        PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.02]}, 
        DisplayFunction -> Identity]

plt7 = 
      ListPlot[Transpose[{y, (Subscript[μ, 2  ex])}], PlotJoined -> True, 
        PlotStyle -> {RGBColor[0, 1, 0], PointSize[0.02]}, 
        DisplayFunction -> Identity]

Show[plt6, plt7, DisplayFunction -> $DisplayFunction]