Hydrocarbon Toolkit in Action

The following describes the hydrocarbon toolkit, its properties, and some simple calculations that can be made using the toolkit

Copyright 2000 Peter Meulbroek

The following assumes some degree of awareness of Mathematica.  While this is probably not the case with the average geochemist, the Mathematica Book is approximately 1400 pages of joy, and can be used to learn simple things.  In order to use the toolkit, you need to be aware of a couple of concepts:
1)  Contexts
2)  Packages

Disclaimer

By using this toolkit, you are implicitly agreeing to the following statement:

1.  Use of the toolkit is automatically granted to individuals.
2.  Changes to source are automatically allowed for personal use, but NOT for distribution.  Send any additions to me, and I will gladly incorporate them.
3.  Distribution of the source is allowed, but only to individuals.  Distribution must be for free;  no fees may be associated with it.
4.  The Toolkit may not be distributed in any other package, or bundled with any other software.
4.  Each file must contain the following statement in a prominant location:
"<filename> was written by Peter Meulbroek, who retains all rights.  Use of <filename> is permitted by individuals.

Setup

The easiest way to use the hydrocarbon toolkit is to implement the following line in your code:

     Map[Needs, {"Phase`Master`", "Standard`Master`", "Graphics`MultipleListPlot`"}] ;

This loads three packages:
    Phase`Master, which defines the hydrocarbon toolkit symbols
    Standard`Master, which defines several useful functions
    Graphics`MultipleListPlot, which is a Mathematica package used for graphing rows of matrices.

General Background

The toolkit is organized into 3 sections. This section talks about the history of the toolkit, the organization of the toolkit

The toolkit was first envisioned in 1992, when Professor Larry Cathles of Cornell University wanted to predict the volume of oil that methane could dissolve as a function of pressure and temperature. Seemed like such a simple question. He set his fresh graduate student Peter Meulbroek, to spend a little time understanding how to predict this stuff. A little time became a lot of time, Peter became an expert in the application of Equations of State to Liquid-Vapor equilibrium in hydrocarbon fluids, (unusual for a geologist, especially in a University without a Petroleum Engineering Dept.), and the problem was solved with the answer "A lot".

The kit is organized to produce 2 types of numbers:
1) Volume/Density of mixtures
2) phase distribution of mixtures

Volume Calculations

Volume calculations are a bit tricky. Due to the requirements of the flash calculations, the toolkit distinguishes between "Liquid Volume" and "Vapor Volume". However, It is not clear a priori if a single phase mixture is best described as a liquid or as a vapor. The absolute correct way to figure it out is to look either 1) where the bubble/dew point is (if the mixture is subcritical) or 2) look at the total gibbs free energy of the fluid in either phase.

First, consider the compressibility of liquid water

Table[p * LiquidVolume[p, t, {1}, {384},
        VanderWaals]/R / t, {p, .001, 20, 1}, {t, 200, 500, 10}] ;

ListPlot3D[%,  ColorFunction -> (Hue[0.9 #] &),

Mesh -> False,

Ticks -> { TakeEvery[NumberList[Table[t, {t, 200, 500, 10}]], 5], Rest[TakeEvery[NumberList[Table[p, {p, 0, 20, 1}]], 5]], Automatic},

PlotLabel -> "Van der Waals Water Compressibility",

AxesLabel -> {"Temperature (摘)", "Pressure (bars)", "Z"}, ViewPoint -> {-2.018, -2.151, 1.659}

[Graphics:images/toolkit_gr_1.gif]

     - SurfaceGraphics -

Next, consider the vapor volume. Note that either volume is not the 'true' volume when the mixture is two-phase (at the saturation line through (p,t) space)

     Table[p * VaporVolume[p, t, {1}, {384}, VanderWaals]/R / t, {p, .001, 20, 1}, {t, 200, 500, 10}] ;

     ListPlot3D[%,  ColorFunction -> (Hue[0.9 #] &),  Mesh -> False, Ticks -> {TakeEvery[NumberList[Table[t, {t, 200, 500, 10}]], 5], Rest[TakeEvery[NumberList[Table[p, {p, 0, 20, 1}]], 5]], Automatic},  PlotLabel -> "Van der Waals Water Compressibility", AxesLabel -> {"Temperature (\272K)", "Pressure (bars)", "Z"}, ViewPoint -> {-2.018, -2.151, 1.659} ]

[Graphics:images/toolkit_gr_2.gif]

     - SurfaceGraphics -

Note that we have rotated the figure for better viewing.

     Show[%, %%%]

[Graphics:images/toolkit_gr_3.gif]

     - Graphics3D -

We can also use another model to generate the same information

     Table[p * LiquidVolume[p, t, {1}, {384}, AasbergPetersen]/(R    t), {p, .001, 20, 1}, {t, 200, 500, 10}] ;

     ListPlot3D[%,  ColorFunction -> (Hue[0.9 #] &),  Mesh -> False, Ticks -> {TakeEvery[NumberList[Table[t, {t, 200, 500, 10}]], 5], Rest[TakeEvery[NumberList[Table[p, {p, 0, 20, 1}]], 5]], Automatic},  PlotLabel -> "A-P Water Compressibility", AxesLabel -> {"Temperature (\272K)", "Pressure (bars)", "Z"} ]

[Graphics:images/toolkit_gr_4.gif]

     - SurfaceGraphics -

Densities are just the molecular weight divided by the specific volume

     Table[MW[[384]]/LiquidVolume[p, t, {1}, {384}, AasbergPetersen], {p, .001, 20, 1}, {t, 200, 500, 10}] ;

     ListPlot3D[%,  ColorFunction -> (Hue[0.9 #] &),  Mesh -> False, Ticks -> {TakeEvery[NumberList[Table[t, {t, 200, 500, 10}]], 5], Rest[TakeEvery[NumberList[Table[p, {p, 0, 20, 1}]], 5]], Automatic},  PlotLabel -> "A-P Water Compressibility", AxesLabel -> {"Temperature (\272K)", "Pressure (bars)", "Z"} ]

[Graphics:images/toolkit_gr_5.gif]

     - SurfaceGraphics -

Note that water is VERY incompressible at STP compared to oils, and so density is pressure insensitive near atmospheric...

We can also do mixtures. Consider a mixture of C1 -C9. We will look at its VAPOR volume.

     idx = Join[{4},    Alkanes[[Range[9]]]] ;

     Speciesnames[[idx]]

     {methane, ethane, propane, n-butane, n-pentane, n-hexane, 
      
        n-heptane, n-octane, n-nonane, n-decane}

     z = Norm[Join[{1}, Exp[-.3 * Range[9]]]]

     {0.2727617891635065`, 0.20206690331807078`, 0.14969484377475792`,
      
         0.11089666781044388`, 0.08215427212686459`, 
      
        0.06086138169842545`, 0.045087220498058425`, 
      
        0.03340143446485591`, 0.02474439124847157`, 
      
        0.01833109589654498`}

     Show[GasChrom[Tb[[idx]], z], Frame -> True, FrameLabel -> {"Tb", "MF"}, PlotLabel -> "Test Oil Composition"]

[Graphics:images/toolkit_gr_6.gif]

     - Graphics -

     Table[(z . MW[[idx]])/VaporVolume[p, t, z, idx, AasbergPetersen], {p, .001, 20, 1}, {t, 200, 500, 10}] ;

     ListPlot3D[%,  ColorFunction -> (Hue[0.9 #] &),  Mesh -> False, Ticks -> {TakeEvery[NumberList[Table[t, {t, 200, 500, 10}]], 5], Rest[TakeEvery[NumberList[Table[p, {p, 0, 20, 1}]], 5]], Automatic},  PlotLabel -> "A-P Oil Density", AxesLabel -> {"Temperature (\272K)", "Pressure (bars)", "\[Rho]"} ]

[Graphics:images/toolkit_gr_7.gif]

     - SurfaceGraphics -

What's this discontinuity in density? The cubic equation only has one root for part of the (p,t) conditions listed above. The 'true' vapor root is shown in the base-level section, while the 'platform' area returns the single root answer. In order to get the 'true' value, we have to get ahead of ourselves, and find the distribution between phases for the mixture.

     << "Phase`EquationofState`AasbergPetersen`"

     Table[{p, t, Flash[p, t, z, idx, AasbergPetersen, MaxIterations -> 1]},    {p, .001, 20, 1}, {t, 200, 500, 10}] ;

Flash returns {{x,y}, V} where x is the normed liquid composition, y is the normed vapor composition, and V is the fraction of the mixture (by mole) in the vapor phase

This next complex statement does the following:
      1)  Determine the appropriate volumes for the x and y defined above
      2)  determine the densities of those phases (density = aMW/volume, aMW = MW species .  composition)
      3)  average the densities according to how much material is present in each phase.
      4)  calculate the compressibility of the TOTAL mixture

     Map[Module[{ v1 = LiquidVolume[#[[1]], #[[2]], #[[3, 1, 1]], idx, AasbergPetersen], v2 = VaporVolume[#[[1]], #[[2]], #[[3, 1, 2]], idx, AasbergPetersen],  mw1 = MW[[idx]] .    #[[3, 1, 1]],  mw2 = MW[[idx]] .    #[[3, 1, 2]],  F = #[[3, 2]] }, ({mw1, mw2} . {1 - F, F})/({v1, v2} . {1 - F, F}) ] &, %, {2}] ;

     ListPlot3D[%,  ColorFunction -> (Hue[0.9 #] &),  Mesh -> False, Ticks -> {TakeEvery[NumberList[Table[t, {t, 200, 500, 10}]], 5], Rest[TakeEvery[NumberList[Table[p, {p, 0, 20, 1}]], 5]], Automatic},  PlotLabel -> "A-P Oil Density",  AxesLabel -> {"Temperature (\272K)", "Pressure (bars)", "\[Rho]"} ]

[Graphics:images/toolkit_gr_8.gif]

     - SurfaceGraphics -

Note that the 'true' density of the multiphase mixture is quite a bit lower than the estimated density using a single phase answer.

P-T Flash Calculations

P-T flash calculations are straightforward (as previewed in the last section).

Consider a Mixture of Methane and 2,2-dimethylPentane. The distribution between liquid and vapor phases can be determined using the "Flash" function. The flash function takes the following options:

Flash[p, t, z, idx, EOS], where p is pressure in bars, t is temperature in ºK, z is the mole fraction array of the mixture, idx is the index of the species, and EOS selects the equation of state that the flash function uses. EOS must be one of {Mathias, PredictiveRKS, AasbergPetersen, or VanderWaals}. Note that the length of z and length of idx must be the same. The function returns a list of three elements: {LiquidComposition, VaporComposition, and vapor fraction}. The compositions are normed vectors of mole fractions, and the vapor fraction is the fraction of the mix (by mole) in the vapor phase.
Flash also takes some options. MaxIterations specifies the number of iterations to use. 10 seems to work well. The convergence routine is "an area of active development".

The following function calls flash using the Mathias EOS, over a temperature range of 150ºK to 550 ºK, and a pressure range from 1 to 181 bars

     dta = Table[Flash[p, t, {.5, .5}, {4, 90}, Mathias, MaxIterations -> 10], {p, 1, 181, 2}, {t, 150, 550, 5}] ;

This plots the vapor fraction of the mix

     Map[Last, dta, {2}] ;

     ListDensityPlot[%,  ColorFunction -> (Hue[0.9 #] &), Mesh -> False, FrameTicks -> {TakeEvery[NumberList[Table[t, {t, 150, 550, 5}]], 10], Rest[TakeEvery[NumberList[Table[p, {p, 0, 180, 2}]], 10]], None, None}, PlotLabel -> "Phase Diagram:  CH4:2,2-dimethylPentane", FrameLabel -> {"Temperature (\272K)", "Pressure (bars)"} ]

[Graphics:images/toolkit_gr_9.gif]

     - DensityGraphics -

The phase envelope looks as expected, but it does show a convergence problem around the critical point (the pertuberence at the upper right of the figure). This error is not unexpected, as critical points represent the toughest convergence. This is "an area of active research".

The same function, but with the Aasberg-Petersen EOS

     dta2 = Table[Flash[p, t, {.5, .5}, {4, 90}, AasbergPetersen, MaxIterations -> 10], {p, 1, 181, 2}, {t, 150, 550, 5}] ;

This plots the vapor fraction of the mix

     Map[Last, dta2, {2}] ;

     ListDensityPlot[%,  ColorFunction -> (Hue[0.9 #] &), Mesh -> False, FrameTicks -> {TakeEvery[NumberList[Table[t, {t, 150, 550, 5}]], 10], Rest[TakeEvery[NumberList[Table[p, {p, 0, 180, 2}]], 10]], None, None}, PlotLabel -> "Phase Diagram:  CH4:2,2-dimethylPentane", FrameLabel -> {"Temperature (\272K)", "Pressure (bars)"} ]

[Graphics:images/toolkit_gr_10.gif]

     - DensityGraphics -

Converted by Mathematica May 6, 2000