Solubility of Methanol in Natural Gas
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)^5) 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)^5) 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)^5)(∑_(j = 1)^5) 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)^5) y[i] Bp[i];

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

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

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

k[1, 1] = 0; k[2, 2] = 0; k[3, 3] = 0; k[4, 4] = 0; 
    k[5, 5] = 0;
k[2, 1] = k[1, 2]; k[3, 1] = k[1, 3]; k[4, 1] = k[1, 4]; 
    k[2, 3] = k[3, 2]; k[3, 4] = k[4, 3]; k[2, 4] = k[4, 2]; 
    k[1, 5] = k[5, 1]; k[2, 5] = k[5, 2]; k[3, 5] = k[5, 3]; 
    k[4, 5] = k[5, 4];
k[1, 2] = -0.035; k[1, 3] = -0.2141; k[1, 4] = 0; 
    k[3, 2] = 0.035999; k[4, 2] = 0.020997; k[4, 3] = 0.0115; k[5, 1] = 0; 
    k[5, 2] = 0.202; k[5, 3] = -0.035999; k[5, 4] = 0.02530;

R = 1.987;

Critical Pressure and Temperature and acentric factor for Methanol, CH4, N2, CO and H2

Pc[1] = 1153.6; Tc[1] = 923.7; ω[1] = 0.5556;

Pc[2] = 673.1; Tc[2] = 343.9; ω[2] = 0.0;

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

Pc[4] = 507.4; Tc[4] = 239.3; ω[4] = 0.048;

Pc[5] = 190.8; Tc[5] = 59.7; ω[5] = 0;

Flash Temperature is varied from  -90 °C to 40 °C and Pressure is 10 bars. 
Feed is arbitrarily composed of 50% Methanol and 50% Natural Gas (70% CH4, 10% N2, 10% CO and 10% H2).

P = 14.696 10;

z[1] := 0.5; z[2] = 0.5 0.7; z[3] = 0.5 0.1; z[4] = 0.5 0.1; 
    z[5] = 0.5 0.1;

{-90, -80, -70, -60, -50, -40, -30, -20, -10, 0, 10, 20, 30, 40} 1.8 + 32 + 459.67

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[1] -> 0.9, x[2] -> 0.1, x[3] -> 0.01, 
        x[4] -> 0.01, x[5] -> 0.01, (Subscript[ϕ, vap]) -> 0.5, 
        y[1] -> 0.00001, y[2] -> 0.65, y[3] -> 0.1, 
        y[4] -> 0.1, y[5] -> 0.1, (Subscript[Z, l]) -> 0.1, 
        (Subscript[Z, v]) -> 0.9}

For[i = 1, 
      i < 15, {T =  {329.67, 347.67, 365.67, 383.67, 401.67, 419.67, 437.67, 
              455.67, 473.67, 491.67, 509.67, 527.67, 545.67, 
              563.67}[[i]], 
        sol[i] = 
          FindRoot[{Eq[1], Eq[2], Eq[3], 
                Eq[4], (∑_(i = 1)^5) y[i] == 1, EQ[1], EQ[2], 
                EQ[3], EQ[4], (∑_(i = 1)^5) x[i] == 
                  1, (∑_(i = 1)^5)(( 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]}, {x[3], 
                x[3] /. sol[i - 1]}, {x[4], x[4] /. sol[i - 1]}, {x[5], 
                x[5] /. sol[i - 1]}, {(Subscript[ϕ, vap]), (Subscript[ϕ, vap]) /. 
                  sol[i - 1]}, {y[1], y[1] /. sol[i - 1]}, {y[2], 
                y[2] /. sol[i - 1]}, {y[3], y[3] /. sol[i - 1]}, {y[4], 
                y[4] /. sol[i - 1]}, {y[5], y[5] /. 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++}]

T = {-90, -80, -70, -60, -50, -40, -30, -20, -10, 0, 10, 20, 30, 40}

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

tbl = Transpose[{T, Y}]

<< Graphics`Graphics`

LogListPlot[tbl, PlotStyle -> RGBColor[1, 0, 0], 
      PlotJoined -> True, PlotRange -> {10^-7, 1}, 
      AspectRatio -> 1.15]

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