version 1.7.0.2 (75.6 KB) by
Mark Mikofski

water and steam properties and derivatives for MATLAB

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'

Mark Mikofski (2021). IAPWS_IF97 (https://github.com/mikofski/IAPWS_IF97), GitHub. Retrieved .

Created with
R2007b

Compatible with any release

**Inspired by:**
X Steam, Thermodynamic properties of water and steam.

**Inspired:**
Water and Steam Refractive Index, A Simple Finite Volume Solver for Matlab

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

jinshah kalluvilaMax FeineisP AStanislav JasoAny chance to also include entropy in the list of functions? This is very important for turbine models

mai emyhow to use it in Matlab?

Mark Mikofski@Christian oops, sorry! I got water and ideal gases confused for a second! See what happens when you start counting photons instead of molecules? So nope, I don't think using `cv=R-cp` would be correct, but maybe you can use `h=u+pv` and then `cv=du/dT` at constant volume. Ideally this should be resolved by a PR against issue #5.

Mark MikofskiHi Christian, Thanks for your comment. The specific heat at constant pressure as a function of pressure & enthalpy, `cp_ph()` is already explicitly available. Is this sufficient for your needs? I just clarified a note in the help string and README, where there were some typos, but the calculation code is unchanged. Try: `cp = IAPWS('cp_ph', 30, 2000)` as an example of `cp` at 30[MPa] and 2000[kJ/kg]. As for specific heat at constant volume, there is an open issue to add this (https://github.com/mikofski/IAPWS_IF97/issues/5), but as a workaround, I think you can use `cv = R - cp`, and I think you can get `cp_pT` by combining `cp_ph` and `h_pT`. Hope that helps. Honestly, since I've switched from solar thermal to PV and from MATLAB to Python, I haven't been maintained this, so if you are interested in taking over this GH repo or this MATLAB FEX, you're welcome. Otherwise, please consider submitting a PR and I will merge it. Cheers!

Christian SoulardHello Mark,

Is there any reason for not having explicit functions for specific heat capacities (cp and cv) ? Agreed, they can be calculated as a function of derivative of other thermodynamic properties but an explicit function would be useful for the user. Xsteam has them implemented (cp and cv at the saturation and as a function of various input combinations).

Thanks for the great work !

zhangfu xuMark MikofskiFYI: I've pushed an update that removes the @folder to a separate branch, so hopefully that solves this one issue. BTW: this also works just fine Octave.

Mark MikofskiJan, Sorry for the inconvenience. Please either delete or rename the folder called "@IAPWS_IF97" it should not be there. I have just cloned the repository, and received a similar error which I traced to the erroneous folder. Renaming the folder from @IAPWS_IF97 -> _IAPWS_IF97 solves the problem. Again my apologies

Jan HeitlandHi Mark, thanks for your work,

i've got an issue with the example file:

Error using IAPWS_IF97

Too many input arguments.

Error in IAPWS_IF97_example (line 10)

h = IAPWS_IF97('h_pT',p,T); % [kJ/kg] enthalpy = f(p,T)

Thanks for your help

samsamJITENDRA KEWALRAMANIHi,

I am very new to Matlab and IAPWS formulation too

I was trying to have IAPWS 2016 Formulation into MatLab. I did basic coding for Helmholtz free energy, having 2 parts free energy and residual energy, as per IAPWS 2016 Version.Link for getting coding(document) is as follow :( https://njit0-my.sharepoint.com/personal/jak93_njit_edu/_layouts/15/guestaccess.aspx?docid=1b9b6fd983fc44bca8474e0102985f525&authkey=ASEeFJ-6pK1fwxNx0S8QBlA&e=uR4N77)

In my coding all the free energy,derivative,constant,required for calculation are coming perfect( I checked on online calculation page of IAPWS-2016) But when it comes to summation of residual energy final output of matlab and IAPWS does not match.

It would be very helpful if someone can help me please.

vaya putrahi

i am using Matlab R2017a

when i try this, some mistake, anyone help me

Not enough input arguments.

Error in XSteam (line 221)

fun=lower(fun);

Peter CookMark,

I did some forking of your package last week - nothing major but I made the following properties available as low-level calls:

'd2gammadtau2_pT','d2gammadpi2_pT','d2gammadtautau2_pT','d2gammadpipi2_pT','d2gammadpitau2_pT',...

'w1_pT','w2_pT','u1_pT','u2_pT','cv1_pT','cv2_pT',...

'gamma1_pT','gamma2_pT','s1_pT','s2_pT','phi3_rhoT'

The updated region 2 dimensionless Gibbs derivatives (d2gammad...) return the ideal and residual derivatives separately - I'd like to discuss if this is actually necessary as IAPWS IF97 Table 13 leads me to believe some of these calls are redundant.

w1_pT & w2_pT return the region 1 & 2 sonic velocities. I have an outside function that computes a retarded 2 phase/equilibrium sonic velocity that uses the same approximation as FLUENT, but that property (2-phase speed of sound) is not supported by the IAPWS - let me know if you are interested in that as well.

u1 & u2 return the region 1 & 2 internal energy.

cv1 & cv2 return the region 1 & 2 isochoric heat capacity.

gamma1 & gamma2 return the region 1 & region 2 dimensionless Gibbs energy.

s1 & s2 return the region 1 & region 2 entropy.

phi3_pT returns the region 3 dimensionless Helmholtz energy. I stopped here because I realized I very rarely use region 3 properties in my own work.

Send me an email and we can discuss a couple of these - I can add this to the Git for the code as well if you'd like.

Maroua ROUABAHGreat job!

Please, can this program be extended for calculating thermodynamic properties of ice VII (from 2 to 20GPa) beyond 1000 bar?

Haben Ghebremedhin@ Mark Mikofski Outstanding job on vectoring the Xsteam file, I was wondering if you have any suggestion on how to calculate Internal energy, Isotherm bulk modulus and coefficient of thermal expansion using your function?

It there is a way to calculate these values using your function, it would be a perfect fit for the Thermal Liquid Settings module in simscape.

Leonardo PaoliHi, Amazing work!

However I am having an issue, with temperatures above 800 C.

For example: IAPWS_IF97('T_ph',20,4300) = NaN

when its supposed to simply give 1164.20 K

I had the same problem with XSteam

Mark MikofskiPeter, good catch! The issue is not actually with k_pT, but with psat_T(Tsat_p(p)) is not exactly equal to p, due to floating point errors.

Unfortunately, not sure of the best fix for it, but probably a tolerance maybe based on eps(p)

I think the culprit is on line 159:

https://github.com/mikofski/IAPWS_IF97/blob/master/IAPWS_IF97.m#L159

valid1 = p>=psat & p<=pmax & T>=Tmin & T<=TB13; % valid range for region 1, include B13 in region 1

since psat is defined as psat_T(T), it should return exactly what was put in for p, ie: linspace(3, 4, 10), but floating point errors are confounding this, and we get::

K>> p>=psat

ans =

false, false, true, true, false, true, false, false, false, false

K>> format long

K>> err = abs(p-psat)

err =

1.0e-13 *

0.297539770599542 (fails)

0.124344978758018 (fails)

0.133226762955019

0.013322676295502

0.093258734068513 (fails)

0.008881784197001

0.031086244689504 (fails)

0.532907051820075 (fails)

0.275335310107039 (fails)

0.266453525910038 (fails)

You can see that the numbers, due to floating point errors are off by less than 1e-12, but eps(p) is ~ 1e-15, so if I use 100*eps(p) it works:

K>> abs(p-psat) < eps(p*100)

returns all true. If you think this is a good fix, submit it as a PR on GitHub, I will merge it add you as author and push it up to MATLAB FEX. Thanks for your contribution!

Mark MikofskiHi Peter, thanks for your comment. This MATLAB IAPWS97 implementation is on GitHub here:

https://github.com/mikofski/IAPWS_IF97.

Can you please submit an issue on GitHub:

https://github.com/mikofski/IAPWS_IF97/issues

Also feel free to hack on the code and send a PR. I will update the MATLAB FEX version with your changes and credit them to you.

Peter CookThis is a magnificent package, though I've got a small quirk to report for 'k_pT' calcs on the 1-2 boundary (region 4) - it appears the function is only correctly identifying the saturation line about 70% of the time. I know this isn't a big deal physically because for most heat transfer problems here most energy is lost through latent heat, but I thought it was worth reporting. Here's an example:

>> P = linspace(3,4,10); k = IAPWS_IF97('k_pT',P,IAPWS_IF97('Tsat_p',P))

k =

0.046698

0.047223

0.63122

0.62921

0.048767

0.62522

0.049777

0.050277

0.050774

0.051268

Peter PauskaKymbat AtamuratovaMark, good day.

I have problem with this program...

When compiling says

Error in ==> IAPWS_IF97_example at 10

h = IAPWS_IF97('h_pT',p,T); % [kJ/kg] enthalpy = f(p,T)

How may it be corrected? "Pressure and temperature ranges" work,but enthalpy is not determined. Thanks!

http://www.mathworks.com/matlabcentral/fileexchange/35710-iapws-if97-functional-form-with-no-slip/content/html/IAPWS_IF97_example.html

NahlaMany thanks Mark

Mark MikofskiSo Sorry @Nahla. You are absolutely correct! I'm sorry I keep inverting Cp by mistake!

It is defined exactly as you have stated it and as it should be according to thermodynamics, the IAPWS and the wikipedia link I originally pasted:

cp = (dh/dT)|p

So, so, so sorry for the confusion!

NahlaOh ok, then probably it is the way you defined it in the program, because the heat capacity definition is

Cp=(dH/dT)|p

Thank you Mark, I appreciate it.

Mark Mikofski@Nahla - yes exactly, I am using the relation:

Cp(p, h) = (dT(p, h)/dh)|p

everywhere in IAPWS_IF97

NahlaThank you again Mark, sorry if my question sounds silly, but do you mean cp_ph is the reciprocal of the derivative of the temperature of steam wrt enthalpy? because this is what it is according to the thermodynamic relations.

Mark Mikofski@Nahla, thank you very much, I'm glad you've found this MATLAB implementation of IAPWS-IF97 useful.

The derivative of the temperature of steam wrt enthalpy is equivalent to cp_ph. Hope that helps!

See http://mikofski.github.io/IAPWS_IF97/ and https://en.wikipedia.org/wiki/Heat_capacity#Thermodynamic_relations_and_definition_of_heat_capacity

NahlaThis is a very useful code that I have been using now for years, many thanks Mark.

I have a question, is there a way I can find the derivative of the temperature of steam wrt enthalpy?

Mark MikofskiCharlie 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?

Mark MikofskiFernando 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

Fernando BelloI 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

Charlie HansenMark,

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.

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

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

Mark Mikofskiok 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!

Mark MikofskiNahla, 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!

NahlaI 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

EhsanZoefast nice... parallel computing :)

Mark MikofskiThe documentation is online here:

http://mikofski.github.io/IAPWS_IF97/

Thomas ClarkNice vectorisation, thank you for this Mark.

Peter SugimuraMark Mikofskiuntil 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)

Mark MikofskiGreat 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.

Jonathan CurrieGood 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.

Mark MikofskiI'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?