(Redirected from Richard Towers)
Jump to: navigation, search

Richard Towers is a GSOC2011 student (Masters of Engineering, Mechanical) from The University of Durham (UK), working on improvements to FPROPS.

Browse code: richard:models/johnpye/fprops


Goals have been updated as of 4 August 2011 and are shown below. The initial project goals are not vastly different, but can be seen in my GSOC application.

  1. Complete integration of PR EOS in FPROPS. All the departure functions should be implemented, tested, working, and connected with the ideal gas functions, for calculation of all the 'raw' properties complete and correct.
  2. PR EOS should be tested with the current sat.c Akasaka approach for calculating single-fluid saturate mixture states. Algorthim should be working and giving reasonable values, and reasonable error messages if failing.
  3. Existing Helmholtz code should be modified to use the new data structures and all existing tests should still be working correctly with no changes having arising from the code refactor.
  4. Existing Python and ASCEND bindings need to be updated to work with all the new code, and be operating correctly.
  5. Existing build tools should be updated to operate on the new code, and FPROPS built and tested using SCons on all three operating systems.
  6. Progress on XML and XSLT tools for FPROPS should be documented, with remaining tasks/intentions detailed.
  7. FPROPS wiki page updated with information about the modifications completed, new functionality, etc.
  8. Convert data from RPP book to C data structures to increase fluid support to 50+ fluids.

Progress reports

17th August 2011

Peng robinson and helmholtz in agreement on pressures now. If I can just get the helmholtz energy departure function working then it'll be working in time for the deadline.

12th August 2011

Difficult week, main computer broken, difficulties compiling on second. Finally succeeded in running scons. Work on python bindings, they're working now but need tidying up.

Thought for the day: SWIG is really clever.

4th August 2011

Helmholtz calculations now work with the new data structures. All fluids and tests have been updated and work (except for ethanol and isohexane who's tests report small and large errors respectively). Peng-Robinson appears to give almost completely different values to helmholtz even for pressure which should be trivial. This is not expected, and the cause needs investigating before we can move into more FPROPS-esque support for ideal and saturation calculations. John emailed with an updated list of goals which are updated above.

25th July 2011

Peng-Robinson EOS now successfully calculates p(V,T), V(p,T) and \Delta h(V,T), \Delta s(V,T). Redlich-Kwong EOS also works, but has additional untested support for a,u,g departure functions. Once the two new EOSs are integrated with current ideal and saturation calculations then we'll have two new EOSs.

15th July 2011

Most of the way to a working implementation of the Peng-Robinson equation of state, there are still some issues to be resolved with respect to ideal properties and internal representation of data.

XML data file structure decided upon and a human readable transformation written, still require an input form to create these files. Will also need to modify cToXML to create files for the existing fluids (Note: Will have to be careful with rounding).

4th July 2011

Having discussed flexibility with regards to equation of state decided on the future structure of the program. Using various structs and unions fprops now has a skeleton capable of supporting several possible equations.

Still to do in this area: Actually implement other EOSs (peng-robinson first), Refactor FPROPS to fit with new way of doing things.

24th June 2011

Switched to Mac OS X for development, switched to the Xcode IDE for editing, hopefully the overheads in learning the new system will be worth it. Spoke to John about support for multiple equations of state, created a couple of files ready to start experimenting.

20th June 2011

All fluid files converted from source to XML, parser demonstrated to work for every file tried so far. Further features (multiple correlations from multiple sources etc.) should be relatively easy to implement now. Acquired a copy of Sandler "Chemical and Engineering Thermodynamics" which will be helpful in sorting out implementations of the cubic equations of state.

17 June 2011

Finished writing an xml parser and translating c data files to xml .fluid files. readFluidFile() from parser.c will now load data at run time. There's still a few features to be implemented here but it's up and running now.

12 June 2011

Found a function to calculate critical pressure in sat.c, made this public and used it in redkw.c. In the process of testing the results now. Should still probably rewrite in terms of reduced vars, and still need to think about V_0. Also remaining dep. functions need deriving. Once that's nailed there's some computing on the horizon...

11 June 2011

The simpler departure functions for Redlich Kwon are nearly implemented now (see redkw.c). The functions should probably be rewritten in terms of reduced variables, which may eliminate the (currently terminal) problem of unspecified critical pressure. Also careful consideration needs to be given to the value of V^0, this seems to be related to a reference state, which needs to be defined. Finally the calculation of Redlich-Kwong coefficients should be moved so that they're calculated a minimum number of times, we'll need some variables for this.

9 June 2011

See Redlich Kwong EOS in FPROPS for details of current thoughts on its implementation.

4 June 2011

Record of an email sent by John Pye to Richard:

Looking at Reid, Prausnitz & Poling (I have the 4th edition) it might be that you will have an easier time getting started with Redlich-Kwong equations of state, rather than the Peng-Robinson ones. The EOS is detailed in Section 3-4 "Cubic Equations of State" and the derivation of so-called 'departure functions' is detailed in Example 5-1 "Derive the departure functions for a pure material...". This derivation includes functions for calculation of extensive properties A, S, H, U, G. Some additional partial derivatives should give you Cp, Cv, and then you need also to convert to the intensive properties. How to use these departure functions to calculate is detailed in Section 5-2 "Fundamental Thermodynamic Principles", eg eq 5-2.4.

Let me put the above into the context of the current code. The eq 5-2.5 contains an ideal part (an integral of Cp°), a 'departure function'. The departure functions are functions exclusively of the p-v-T surface, ie the EOS. So given the EOS we can do the maths to work out what these departure functions need to be. This is done for the Helmholtz approach in the file helmholtz.c: see the functions eg helmholtz_h_raw, which shows the ideal and the 'residual' (departure) parts clearly. The residual parts are then implemented further down in that file.

I would propose then that to do a RK or PR EOS, we can keep most of helmholtz.[ch] and simply calculate new functions like helm_resid_tau and so on, specific to the EOS in question. All the current stuff from ideal.c continues to be used as-is, and the saturation calculation also remains as-is. We just need new functions like redkw_resid, redkw_resid_del, redkw_resid_tau, etc.

Doing the maths for those functions will need to be extremely carefully done, as from my experience is very very easy to introduce little calculation errors along the way. A good idea is to make use of tools like wxMaxima or similar, to do some of the maths symbolically, even if only as a double-check for hand-calculation. You then need to take the resulting expressions and work out how to code them efficiently.

We will of course then need to work out suitable data-structrues that allow functions like helmholtz_h_raw (renamed to fprops_h_raw, perhaps) to be common. That will require the HelmholtzData structure to be generalised to contain some pointer to the a struct containing the necessary residual functions, I suspect.

Further to the above, we will have to ponder how to implement support for thermo derivatives (derivs.[ch]). I think that this requires additional functions helmholtz_betap and helmholtz_alphap to be generalised as well. That might be trivial, actually.


Need to find correct data for comparison. Working version of REFPROP? Literature? Inconsistencies need to be investigated and hopefully fixed.

Improvements to documentation

As part of the build up to starting work I've been reading through the LyX generated help file that comes with ASCEND, Having just finished a dissertation written in LyX and LaTeX/pgf/tikz I could probably help significantly with the general appearance of this document, for example:

Error creating thumbnail: File missing
Old Types Diagram
Error creating thumbnail: File missing
New Types Diagram
An old diagram showing a thin walled cylinder
Old Cylinder Diagram
A new diagram showing a thin walled cylinder
New Cylinder Diagram

Separating Correlation Coefficients

Coefficients should be stored in human readable files, not as C source. The plan is to create an XML structure for these files and a function to parse the files, loading fluid information at runtime. John suggests using RelaxNG as the schema, I'm not sure whether or not ASCEND already has a preferred cross platform XML parser, to begin with I'll use GNOME's libxml.

Relax NG Grammar

Element names identify the property in question and attributes identify values. The advantage of using attributes instead of node text is that we can specify different attributes (for the same property) that are treated differently, for example: it is sometimes useful to specify the constant ideal term as a multiple of the specific gas constant instead of as an absolute value. fluidSchema.rng

Example XML Fluid File

As an example a fluid file for Hydrogen can be found in hydrogen.fluid. Note that the constant term is specified as a multiple of R as described above.


The minimum requirements of the parser are quite simple, however it would be nice not to have to worry about arranging xml elements in a specific order and some variables are easier to specify as multiples of some other variable. The beginnings of a competent parser are found in parser.c.

Equations of State

Firstly I'll be looking at the Redlich Kwong & Peng Robinson cubic EOSs. See richard:models/johnpye/fprops/pengrob.c, richard:models/johnpye/fprops/pengrob.h, richard:models/johnpye/fprops/redkw.c and richard:models/johnpye/fprops/redkw.h for more. Later an attempt might be made at implementing Lee and Kesler's modified Benedict-Webb-Rubin EOS, but this is another job for another day...

One of the main challenges here is to write functions which can take any correlation data struct as an input, decide which correlation to use, and return the desired value. After some consultation with John, we've decided on a struct containing all of the data and pointers to the various functions. A union is used to enable different types of data to be passed to functions with the same arguments.

/** Enumerate EOS correlations
enum eos{
    mbwr //etc.

/** Union for equation of state data
typedef union eos_tag{
    HelmholtzData* helmholtzData;
    CubicData* cubicData;
} EosData;

/** Definition of a fluid property function pointer
typedef double (*propFnPtr)(double,double,const EosData*);

	Structure pointing to the necessary data and functions to calculate
    fluid properties.
typedef struct pureFluid{
	enum eos type;
    IdealData *idealData;
	EosData  *eosData;
	//Pointers to departure functions
    propFnPtr p_fn;
    propFnPtr u_fn;
    propFnPtr h_fn;
    propFnPtr s_fn;
    propFnPtr a_fn;
    propFnPtr cv_fn;
    propFnPtr cp_fn;
    propFnPtr w_fn;
    propFnPtr g_fn;
} pureFluid;

The structure can be created as follows, and then the function pointers can simply be followed to the correct function.

EosData eosData={.cubicData=&cData};
    pureFluid fluidData={
        //Pointers to departure functions

double p=fluidData.p_fn(T,rho,fluidData.eosData);