Code covered by the BSD License  

Highlights from
IAPWS_IF97 functional form with no slip

4.875
4.9 | 8 ratings Rate this file 44 Downloads (last 30 days) File Size: 78.9 KB File ID: #35710 Version: 1.7
image thumbnail

IAPWS_IF97 functional form with no slip

by

Mark Mikofski (view profile)

 

18 Mar 2012 (Updated )

Water and steam properties and derivatives based on the IAPWS IF97. Functional form. No slip.

| Watch this File

File Information
Description

IAPWS_IF97(FUN,IN1,IN2) is 27 functions of water properties and derivatives, based on the International Association on Properties of Water and Steam (http://www.iapws.org). Thermodynamic, hydrodynamic and non-linear modeling often requires thermodynamic derivatives, therefore IAPWS_IF97 can calculate most property derivates as functions of pressure and enthalpy, e.g.: dT/dp_ph, cp_ph, dv/dp_ph and dv/dh_ph. Since modeling often involves multiple dimensions that are discretized or meshed to form a set of either finite-difference or finite-element equations, IAPWS_IF97 is vectorized even across regions (subcooled/compressed-liquid, saturated, superheated and supercritical). For arrays is an order of magnitude faster than XSteam and is multithreaded if your computer is capable.
This is the functional form. I will submit a class & package versions definition soon, that also offer slip correction using Zivi's correlation (1964) for 2-phase flow.

Please report bugs here or at Github:
https://github.com/mikofski/IAPWS_IF97/issues

Reference Documents:
Industrial Formulation 1997 (IF97-Rev, IAPWS-IF97), IAPWS-IF97-S01, IAPWS-IF97-S03rev, IAPWS-IF97-S04, IAPWS-IF97-S05, Revised Advisory Note No. 3 Thermodynamic Derivatives from IAPWS Formulations 2008, Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance, 2008 Revised Release on the IAPWS Formulation 1985 for the Thermal Conductivity of Ordinary Water Substance.

Functions:
Quality:
    'x_ph','x_hT','x_pv','x_vT'
Thermal Conductivity [W/m/K]:
    'k_pT', 'k_ph'
Viscosity [Pa*s]:
    'mu_pT', 'mu_ph', 'dmudh_ph', 'dmudp_ph'
Enthalpy [kJ/kg]:
    'h_pT', 'hL_p', 'hV_p', 'dhLdp_p', 'dhVdp_p'
Specific Volume [m^3/kg]:
    'v_pT', 'v_ph', 'vL_p', 'vV_p', 'dvLdp_p', 'dvVdp_p', 'dvdp_ph', 'dvdh_ph'
Temperature [K]:
    'T_ph', 'dTdp_ph', 'cp_ph'
Saturation Pressure [MPa] and Temperature [K]:
    'psat_T', 'Tsat_p', 'dTsatdpsat_p'

Acknowledgements

X Steam, Thermodynamic Properties Of Water And Steam. inspired this file.

This file inspired Water And Steam Refractive Index and A Simple Finite Volume Solver For Matlab.

Required Products MATLAB
MATLAB release MATLAB 7.5 (R2007b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (18)
25 Mar 2015 Mark Mikofski

Mark Mikofski (view profile)

Charlie Hansen: sorry for not seeing your message here. I'm sorry you are running into an issue, but it will be difficult for me to debug it without more info. Can you please post an issue on GitHub here: https://github.com/mikofski/iapws_if97/issues describing in detail what exact code you have tried, what the exact stack trace was use lasterror MATLAB builtin and exactly how you have IAPWS_IF97 on your path?

Comment only
24 Mar 2015 Mark Mikofski

Mark Mikofski (view profile)

Fernando Bello: In order to use `pT_uv()`, which is a reverse function, you will need my `newtonraphson()` solver from the FEX: http://www.mathworks.com/matlabcentral/fileexchange/43097-newton-raphson-solver. Once downloaded either add it to your path or extract to your main MATLAB working folder.
AFAIK `pT_uv()` is the only function that depends on `newtonraphson()`. Note it is *not* an official IAPWS-IF97 function, although it does use the IAPWS-IF97 correlations.
Glad you are finding it useful. Please don't hesitate to report any other issues! --Mark

Comment only
24 Mar 2015 Fernando Bello

I have doubts with this program... When you use

pT_uv(41990.5,1/999.794) it says this error:

Undefined function 'newtonraphson' for input arguments of type
'function_handle'.

Error in pT_uv (line 66)
[x,resnorm,F,exitflag,output,jacob] = newtonraphson(fun,x0,options);

What is it? do i need the newton-raphson solver that is not included in this program.

About the rest of the program, it works great! it is a really good and fine job! Thank you for sharing

08 Jan 2015 Charlie Hansen

Mark,

This seems to be an extremely useful functionality (especially because it is not self-referencing), but I am unable to run it in even a simple test. Trying to run the example file (IAPWS_IF97_example.m) and making my own simple call to the function both result in the same error:

Error using IAPWS_IF97
Too many input arguments.

Is there some incompatibility that you know of that may be causing this issue? I did not alter the main IAPWS_IF97.m file and also tried to download it again but saw no change.

Thanks for your help.

Comment only
20 Sep 2014 Nahla

Nahla (view profile)

Thank you Mark, it is working now. I appreciate it.

Comment only
20 Sep 2014 hgghh

hgghh (view profile)

https://www.drupal.org/u/hd.-clemson.-vs.-florida.-state.-live.-stream.

19 Sep 2014 Mark Mikofski

Mark Mikofski (view profile)

ok Nahla, I found the bugs and fixed them enough so that you should get the correct derivatives for specific volume along the saturated liquid and vapor lines, but you will need to use the new dvLdp_p and dvVdp_p functions. There are 2 new tests included which will plot the derivatives together with finite difference approximations. For more details please see IAPWS "Revised Advisory Note #3, Thermodynamic Derivatives from IAPWS Formulations" (http://www.iapws.org/relguide/Advise3.pdf) which details how derivatives are found. To get the derivatives on the liquid and vapor lines, since the specific volume and enthalpy are functions of p or T only then a Taylor expansion, dhL/dp = (dh/dp)_T + (dh/dT)_p * dT_sat/dp_sat is used where the partial derivatives (dh/dp)_T and (dh/dT)_p are found using the advisory note separately for regions 1, 2 and 3. Ditto for specific volume. I made some silly sign errors, and I also switched cv accidentally for cv+p*v*alphap - I think my eyes wandered to internal energy, u, in the table instead of enthalpy, h, which is right below it. Sorry for the inconvenience! Thanks for noticing the error!

Comment only
19 Sep 2014 Mark Mikofski

Mark Mikofski (view profile)

Nahla, I think you may have discovered a typo in IAPWS_IF97. If I change the signs of the region 4b terms in dhLdp on line 914 to ...

dhLdp(valid4b) = (v3L - Tsat4b.*alphap3L./betap3_rhoT(rho3L,Tsat4b).*(1 - p4b.*alphap3L.*dTsatdpsat4b))/conversion_factor + cv3_rhoT(rho3L,Tsat4b).*dTsatdpsat4b;

I get a curve that matches a finite difference approximation of dhLdp. This function is also used in dvdp_ph(). Can you please post an issue on Github, and try to correct the error and validate it versus the IAPWS_IF97 documentation. Thanks!

Comment only
19 Sep 2014 Nahla

Nahla (view profile)

I have a question regarding the derivatives. I need the derivative of density wrt saturation pressure. However when I use the derivative of sp. volume wrt to p and multiply it with -rho^2, it does not give me accurate answers. Is there another way around this?
Thanks

19 Jun 2014 Ehsan

Ehsan (view profile)

 
17 Oct 2013 Zoe

Zoe (view profile)

fast nice... parallel computing :)

14 Oct 2013 Mark Mikofski

Mark Mikofski (view profile)

The documentation is online here:
http://mikofski.github.io/IAPWS_IF97/

Comment only
13 Dec 2012 Thomas Clark

Nice vectorisation, thank you for this Mark.

15 May 2012 Peter Sugimura  
26 Mar 2012 Mark Mikofski

Mark Mikofski (view profile)

until implemented, quality can be calculated from either enthalpy or specific volume and from either psat or Tsat as follows:

x_ph = @(p,h)(h - hL_p(p))./(hV_p(p) - hL_p(p))

then use x = x_ph(p,h)

or:

x_hT = @(h,T)(h - hL_p(psat_T(T)))./(hV_p(psat_T(T)) - hL_p(psat_T(T)))

then use x = x_hT(h,T)

and the very similar for specific volume:

x_pv = @(p,v)(v - vL_p(p))./(vV_p(p) - vL_p(p))

then use x = x_pv(p,v)

or:

x_vT = @(v,T)(v - vL_p(psat_T(T)))./(vV_p(psat_T(T)) - vL_p(psat_T(T)))

then use x = x_vT(h,T)

Comment only
26 Mar 2012 Mark Mikofski

Mark Mikofski (view profile)

Great suggestions Jonathan. I have added them to the issues list on Github here:
https://github.com/mikofski/iapws_if97/issues
Also, I cp is available as a function of p & h (but not p & T as you pointed out), if that would suit your needs.

Comment only
26 Mar 2012 Jonathan Currie

Good implementation and the vectorization works really well. However would suggest adding:

- Entropy calculations
- Cp, Cv calculations as a function of P,T (not constrained by region)
- Inverse calculations for P
- Quality calculations
- More property pairings (e.g. P,T, P,S, P,H, H,S)
- Customizable units?

Not sure on the implementation question. I've always found classes too slow, there seems to be a bit of getting started overhead in every method call.

22 Mar 2012 Mark Mikofski

Mark Mikofski (view profile)

I've thought of a couple of other ways to present this code:
1. as a "+folder package", each function in it's own file; use namespace or import to use functions; e.g. IAPWS_IF.h_pT(0.1,298)
2. a @class folder, consolidate constants as hidden constant properties, functions as methods or static methods, method to set slip if slip is desired for 2-phase flow
3. "+folder package" for class, same as #2 but in a package for import.
Each has pros and cons. What works best for users?

Comment only
Updates
19 Mar 2012 1.2

revise/format description, add more tags

23 Mar 2012 1.3

add link to github

23 Mar 2012 1.4

show units, categorize and format functions

26 Mar 2012 1.5

acknowledge Magnus Holmgren's XSteam

14 Oct 2013 1.6

add quality, update help, add assertions

22 Sep 2014 1.7

fix bugs in dhL/dp and dhV/dp, and add dvL/dp and dvV/dp functions

Contact us