thorade/helmholtzmedia

Modelica library for the calculation of fluid properties from a Helmholtz energy equation of state (EoS).

helmholtzmedia's Issues

inlining and functions

Some of the functions could be inlined if the intermediate protected variables were used directly; makes the code less readable, but possibly faster.
Some of the more complicated intermediate variables could be calculated in a function that again can be inlined.

upgrade to MSL4

it has been released long time ago

Dependency on Modelica.Media.Examples.Tests

Modelica.Media.Examples.Tests got obsoleted in MSL v3.2.3 and will be removed in MSl v4.0.0. They are still availabe in ModelicaTest though.

The difference between HelmholtzMedia and Modelica.Media?

I am planning to choose a base media library for my model, I am interested in this library, but I am not sure what the difference is between HelmholtzMedia and Modelica.Media, could you give me some clue about it?
I read the paper about this library but didn't get a clue. It would be more user-friendly with an illustration about why we need this new library other than Modelica.Media.

Tidy up HTML in documentation

See https://build.openmodelica.org/Documentation/tidy.filtered

Bridgman 2nd order

A couple of second order derivatives are calculated already, additionally it would be good to calculate the derivatives used by Bridgman for his 2nd order derivatives transformation table. These are:

• (∂²v/∂p²)_T
• (∂²v/∂p∂T)
• (∂²v/∂T²)_p
• (∂c_p/∂T)_p = (∂²h/∂T²)_p

Other base sets would be possible (and easier to calculate from an Helmholtz EoS), but these are the ones suggested by Bridgman.
http://dx.doi.org/10.1103/PhysRev.3.273

Use of HemholtzMedia in openmodelica 3.2.3

Hi, I have downloaded the latest version from GITHUB and loaded the library into OMEdit using 3.2.3. However, none of the examples for tests run. I think I am missing something. If I use ButaneTestModel_ph I get

Check of HelmholtzMedia.Examples.MediaTestModels.ButaneTestModel_ph completed successfully.

Instantiation of HelmholtzMedia.Examples.MediaTestModels.ButaneTestModel_ph completed successfully.

But running fails, the first error is :

Internal error No support of solving not real variables with a non-linear solver.

Equation:
HelmholtzMedia.Examples.MediaTestModels.ButaneTestModel_ph.Medium.SaturationProperties(fixedMassFlowRate.medium.sat.psat, fixedMassFlowRate.medium.sat.Tsat, HelmholtzMedia.Examples.MediaTestModels.ButaneTestModel_ph.Medium.ThermodynamicState(fixedMassFlowRate.medium.sat.liq.phase, fixedMassFlowRate.medium.sat.liq.T, fixedMassFlowRate.medium.sat.liq.p, fixedMassFlowRate.medium.sat.liq.d, fixedMassFlowRate.medium.sat.liq.u, fixedMassFlowRate.medium.sat.liq.h, fixedMassFlowRate.medium.sat.liq.s), HelmholtzMedia.Examples.MediaTestModels.ButaneTestModel_ph.Medium.ThermodynamicState(fixedMassFlowRate.medium.sat.vap.phase, fixedMassFlowRate.medium.sat.vap.T, fixedMassFlowRate.medium.sat.vap.p, fixedMassFlowRate.medium.sat.vap.d, fixedMassFlowRate.medium.sat.vap.u, fixedMassFlowRate.medium.sat.vap.h, fixedMassFlowRate.medium.sat.vap.s)) = HelmholtzMedia.Examples.MediaTestModels.ButaneTestModel_ph.Medium.setSat_p(volume.medium.p) solve for {fixedMassFlowRate.medium.sat.psat,fixedMassFlowRate.medium.sat.liq.h,fixedMassFlowRate.medium.sat.vap.d,fixedMassFlowRate.medium.sat.vap.u,fixedMassFlowRate.medium.sat.vap.h,fixedMassFlowRate.medium.sat.vap.s,fixedMassFlowRate.medium.sat.liq.u,fixedMassFlowRate.medium.sat.vap.phase,fixedMassFlowRate.medium.sat.liq.s,fixedMassFlowRate.medium.sat.vap.T,fixedMassFlowRate.medium.sat.liq.phase,fixedMassFlowRate.medium.sat.liq.T,fixedMassFlowRate.medium.sat.Tsat,fixedMassFlowRate.medium.sat.vap.p,fixedMassFlowRate.medium.sat.liq.p,fixedMassFlowRate.medium.sat.liq.d}

Followed by

Internal error createEquationsForSystems failed then many more

Can you help get me going. Thanks

testing, validation, documentation

issue for collecting all commits that are related testing, validation and documentation

automatic differentiation

• https://github.com/usnistgov/teqp
• https://modelica.org/events/Conference2005/online_proceedings/Session1/Session1c2.pdf
• https://stackoverflow.com/a/64262535/874701
• https://stackoverflow.com/q/78236231/874701
• https://github.com/thorade/HelmholtzMedia/blob/master/HelmholtzMedia/Interfaces/PartialHelmholtzMedium/EoS/setHelmholtzDerivsSecond.mo

TwoPhase_setState fails on Openmodelica 1.18

We've encountered a bug that is also present in HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState. See also #39, but I didn't want to pollute that one with details.

TwoPhase_setState fails in Openmodelica 1.18 (NF) on Windows 10 using either current HelmHoltzMedia master (18ff552) and MSL 3.2.3 or current msl4update (20bd884) and MSL 4.0.0

The error also happens on the OM 1.18 test server (also currently on the OM master):

Regular simulation: ./HelmholtzMedia_HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState -abortSlowSimulation -alarm=480 -lv LOG_STATS
LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
assert            | debug   | HelmholtzMedia_HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState_functions.c:11145: Invalid root: (nan)^(3)
stdout            | warning | Integrator attempt to handle a problem with a called assert.
assert            | debug   | HelmholtzMedia_HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState_functions.c:11145: Invalid root: (nan)^(3)
stdout            | warning | Integrator attempt to handle a problem with a called assert.
assert            | debug   | HelmholtzMedia_HelmholtzMedia.Examples.ConvergenceTest.TwoPhase_setState_functions.c:11145: Invalid root: (nan)^(3)
stdout            | info    | model terminate | Simulation terminated by an assert at time: 10.0002

The test regressed at this point (also cited by @thorade in #39), and I suspect that this OM commit by @perost could be connected/the cause (old trac issue related to that).

The Modelica Specification issue modelica/ModelicaSpecification#2757 is probably connected but does not seem close to resolution.

Now, why all this exposition?
The error seems to occur, as far as we can tell, in HelmholtzMedia.Interfaces.PartialHelmholtzMedium.setState_pd, because

• here f.delta is NaN, with
• delta becoming NaN here during record constructions because it seems to see a d_crit of 0.
• d_crit maybe ends up being 0 a little above when
• Maybe MM is 0 here.

I'm not sure if those 0s occur because something is off in the ordering of binding equations when the function is evaluated, so a default/uninitialized value of 0.0 is used (hence the suspicion of OpenModelica/OpenModelica@14d191b423 above), or because something wonky happens when looking up fluidConstants[1].molarMass.

In any case, the minimal workaround that makes this test (and our internal model) work is to delay the assignment of f.delta:

protected
  //...
  EoS.HelmholtzDerivs f(d=d);
  //...
algorithm
  f.delta := d/d_crit;

This strengthens my suspicion that the ordering of operations could be off in OM.

Questions/remarks:

• Considering the spec issue above, maybe the desired behaviour of this corner of Modelica is not yet fully clear, maybe we should work around/avoid it in the meantime?
• Should I rather report this to OM to fix on their side (e.g. a tweak of OpenModelica/OpenModelica@14d191b423)?
• I'm not competent to assess the wider implications of this workaround, and if this has to be done elsewhere in the library, too. can anybody else contribute?

allow reducing parameters to be different from critical parameters

Necessary for e.g. HMDS Hexa-Methyl-Di-Siloxane.
Some models already support different reducing parameters (e.g. viscosity, vapor pressure, thermal conductivity),
some models need to be changed (e.g. EoS, ancillary density, ...).

HMDS currently uses the reducing parameters everywhere (T_crit is set to T_reduce, even though they are not equal).

setState_dh or setState_pd

Laughman, Qiao, 2016: density as additional state in discrete volumes
https://doi.org/10.1080/13873954.2017.1298625

setState_ph is a 2D iteration. If p,h, and d are known as state, two 1D alternatives would be possible:

• setState_dh
• setState_pd

gas constant

Distinguish between universal gas constant and specific gas constant:
modelica/ModelicaStandardLibrary#2438

Use the identical value for universal gas constant everywhere, instead of best known value at the time the EoS paper was written:
usnistgov/REFPROP-issues#239

Use of uninitialized variable in expression

In this line

HelmholtzMedia/HelmholtzMedia/Interfaces/PartialHelmholtzMedium/Transport/dynamicViscosity_residual.mo

Line 113 in 7b867a5

eta_r := eta_r + visci;

partial values are accumulated into eta_r while executing a for loop. Unfortunately eta_r is uninitialized, so according to Section 12.4.4 this model seems to be invalid.

I would recommend to add an initialization statement eta_r :=0 before entering the for loop.

Non-analytic critical terms for residual Helmholtz energy not yet implemented

Non-analytic terms not yet implemented,
these are needed for the residual Helmholtz energy of water and CO2 in the critical region.

See Span 2000, page 33
http://books.google.de/books?id=dd410GGw8wUC&pg=PA33

Support of HelmholzMedia by OpenModelica

@thorade, we've made good progress with the support of your library in OpenModelica. @perost did an excellent job with the new frontend, so we passed from two models only running with the old front end, to all except one of the 42 runnable models of the library getting successfully compiled into simulation code, see report. Unfortunately, about a quarter of them still fail during simulation because of numerical issues.

I haven't had the time so far to look into that, so I'm not sure whether this is because of some shortcomings of the OpenModelica backend/runtime. Could you have a quick look and report here?

debugging of Ancillary.density_pT_Soave

For Pentane at p=1e-6 and T=578.275 the Soave algorithm goes to the supercritical branch, where three roots exist. According to Span, two of the three roots result in negative densities, but in this case this is not true.
Not a big problem, as this is used as starting value only, but it slows down the pT iteration.

Also, the RKS EoS should be able to correctly describe the liquid region, at least qualitatively.

Ancillary.saturationPressure_T is fixed to model PS5

Currently, only the PS5 model is supported for the ancillary function for saturation pressure. Other models are possible, especially PS6 for water.

setState_hs function not implemented

Adding setState_hs would be helpful for turbine calculation,
when power is given and p_out is to be determined

PIP

G. Venkatarathnama and L.R. Oellrich 2010
http://dx.doi.org/10.1016/j.fluid.2010.12.001

Grüneisen

Grüneisen is used e.g. here:
http://dx.doi.org/10.1063/1.3664084

ɣ = v · (∂p/∂u)_v
ɣ = v/T · (∂T/∂v)_s

time derivatives in two-phase region

HelmholtzMedia/HelmholtzMedia/Interfaces/PartialHelmholtzMedium/specificEntropy_ph_der.mo

Lines 13 to 14 in 258c87e

	s_der := p_der*(-1.0/(state.d*state.T))
	+ h_der*(1.0/state.T);

The time derivative is calculated form the partial derivatives. Here, the single-phase definition was used, but that is not correct for two-phase.

MultiPhaseMixtureMedia working with new interface error

I managed to get MultiPhaseMixtureMedia working but only with the old interface. If I run with the new interface I get the following error

Internal error NFCeval.evalRelationEqual failed to evaluate ‘i == 1‘

Is there any assistance in determining what needs to be fixed to get rid of this type of error in the new interface.

Index reduction does not work as expected

There are some numerical Jacobians when running the Examples.MediaTestModels
with Dymola flag Hidden.PrintFailureToDifferentiate=true
http://www.claytex.com/how-can-i-make-my-models-run-faster/

Where does this have to be fixed?
In BaseProperties, or in density_ph?

top-level folder has to be renamed after downloading & unzipping

After downloading the files and unzipping them, the top-level folder has to be renamed from HelmholtzMedia-master or HelmholtzMedia-0.9.1to HelmholtzMedia, because Dymola (and possibly other tools) expect the top-level folder name to be identical to the package name as specified in the package.mo file.

OpenModelica compatibility

Currently, trying to use HelmholtzMedia in OpenModelica results in errors similar to
Variable fluidConstants[1].molarMass ... and Illegal subscript [1] for dimensions in component ...

fugacity

Add a function that calculates the fugacity
http://books.google.de/books?id=1rOA5I6kQ7gC&&pg=PA243
http://books.google.de/books?id=dd410GGw8wUC&pg=PA40

more fluids

Possibly interesting refrigerants:

• R134a
• R1233zd
• R410A
• R1234ye
• R744
• R32
• R515b
• R1234ze
• R513
• R454b

Possibly interesting reference fluids:

• Argon
• Methane
• SF6
• CO2
• R-134a
• R124
• R125
• Methanol
• Water

Add a new fluid

Hello,

I want to add a new helmholtz media like hydrogen.
I have find these properties in google and I want create a new media.
https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=832374

I don't understand how could I do this ?
Can you help me ?

Thank you

Error using Ammonia Properties

extrapolation

What happens with the EoS and the HelmholtzDerivs at

• T<T_trip
• T>T_max
• d>d_max
• d<0
• p>p_melt ?

See http://dx.doi.org/10.1007/BF02575343
Also see http://dx.doi.org/10.1039/9781849730983-00394
or http://scitation.aip.org/content/aip/journal/jpcrd/44/1/10.1063/1.4913493#s7
http://books.google.de/books?id=-W53DF10lmsC&pg=PA394
http://books.google.de/books?id=-W53DF10lmsC&pg=PA418

How does p_satextrapolate for T>T_crit?
This is also discussed in the dissertation of HD Baehr

cubic EoS

can be transformed into Helmholtz form
http://dx.doi.org/10.6028/jres.121.011

convergence of setSat functions very close to critical point

Very close to the critical point (e.g. T_crit-0.1 or T_crit-1e-6) the setSat functions may fail. This will happen only for fluids where only few data were available for fitting the EoS and the ancillary equations, one example is HDMS.

HelmholtzCool

@CarlosTrujilloGonzalez, awesome work with you fork!
https://github.com/CarlosTrujilloGonzalez/HelmholtzCool

I started something similar here:
https://github.com/UdK-VPT/CoTeTo/tree/master/Generators/Example50JsonFluid
but never finished it. I also started some related discussion here:
CoolProp/CoolProp#115
Happy to discuss how to proceed.

Virial EoS

For the gas region (not liquid!!), a virial equation could be used to calculate start values.
See http://en.wikipedia.org/wiki/Compressibility_factor for an introduction

Add four new properties:
Compressibility factor Z
Second virial coefficient B
Third virial coefficient C
Fourth virial coefficient D

Usually, the input properties are temperature and volume.
Kabelac 1991 shows how to solve for volume, when temperature and pressure are given:
http://dx.doi.org/10.1016/0140-7007(91)90006-3 (in the appendix)

Remaining issues with HelmholtzMedia in OpenModelica

newInst Library report, Library coverage, last 7 runs

---

This is a follow-up ticket for #34, to track the few issues that are left after the release of OpenModelica 1.16.
There is a corresponding ticket in OM Trac: https://trac.openmodelica.org/OpenModelica/ticket/6088

---

1. Examples.ConvergenceTest.Ancillary_Saturation
   This goes into a branch of an if-else that I thought is unreachable (this should also not happen)??? After updating the assert message, it is seen that h=nan

2. Examples.ConvergenceTest.setSat
   An iteration written by me does not converge, close to the critical point, but it works in Dymola. There is a warning about abs value being negative.

3. Examples.MediaTestModels.IsopentaneTestModel (during initialization)
   Failing assert

4. Examples.Validation.Derivatives_SaturationBoundary (during initialization)
   fails during initialization, with liq.s and vap.s having same value

5. Examples.Validation.idealGasLimit
   this should converge to a certain value, but starts to oscillate because of numerical noise, also with other tools, but at a later point. Might be solved by changing some tolerance or other setting?

Helium: transport properties not implemented

Helium does currently use transport properties of Butane!
Instead, it should use the special helium correlations, but these are not yet implemented.

Linux compatibility

Is this library compatible with Linux? I've unsuccessfully attempted to run HelmholtzMedia examples using OMEdit 1.17.0 MSL 3.2.3 on my Linux virtual machine. I've verified that it runs on Windows, but am unsure what the difference is that is causing issues. Thanks.

Example error window:
[1] 09:38:39 Translation Error
[HelmholtzMedia.Interfaces.PartialHelmholtzMedium.Ancillary.saturationPressure_T: 13:3-13:80]: Variable ancillaryCoefficients.pressureSaturation not found in scope Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.Ancillary.saturationPressure_T.

[2] 09:38:39 Translation Error
[HelmholtzMedia.Interfaces.PartialHelmholtzMedium: 571:9-571:52]: Variable sat.liq.d not found in scope Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.setState_pTX.

[3] 09:38:39 Translation Error
[HelmholtzMedia.Interfaces.PartialHelmholtzMedium.EoS.f_i: 9:3-9:64]: Variable helmholtzCoefficients.idealLog not found in scope Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.EoS.f_i.

[4] 09:38:39 Translation Error
[HelmholtzMedia.Interfaces.PartialHelmholtzMedium.EoS.setHelmholtzDerivsSecond: 15:3-15:41]: Class f_i not found in scope Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.EoS.setHelmholtzDerivsSecond (looking for a function or record).

[5] 09:38:39 Translation Error
[HelmholtzMedia.Interfaces.PartialHelmholtzMedium.density_ph_state: 2:1-17:21]: Function .Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.density_ph_der not found in scope Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.

[6] 09:38:40 Translation Error
[HelmholtzMedia.Interfaces.PartialHelmholtzMedium: 732:5-732:62]: Variable helmholtzCoefficients.useLineSearch not found in scope Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.setState_phX.

[7] 09:38:40 Translation Error
[HelmholtzMedia.Interfaces.PartialHelmholtzMedium: 1637:5-1637:96]: Variable sat.vap.s not found in scope Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.saturationTemperature_derp_sat.

[8] 09:38:40 Translation Error
[HelmholtzMedia.Interfaces.PartialHelmholtzMedium.temperature_ph_state: 2:1-18:25]: Function .Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.temperature_ph_der not found in scope Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.

[9] 09:38:40 Translation Error
[HelmholtzMedia.Interfaces.PartialHelmholtzMedium.Ancillary.dewDensity_T: 9:3-9:78]: Variable ancillaryCoefficients.densityVaporModel not found in scope Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.Ancillary.dewDensity_T.

[10] 09:38:40 Translation Error
[HelmholtzMedia.Interfaces.PartialHelmholtzMedium.Ancillary.bubbleDensity_T: 9:3-9:81]: Variable ancillaryCoefficients.densityLiquidModel not found in scope Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.Ancillary.bubbleDensity_T.

[11] 09:38:40 Translation Error
[HelmholtzMedia.Interfaces.PartialHelmholtzMedium.Ancillary.saturationTemperature_p: 16:3-16:80]: Variable ancillaryCoefficients.pressureSaturation not found in scope Modelica.Fluid.Sources.Boundary_pT$boundary1.Medium.Ancillary.saturationTemperature_p.

[12] 09:38:40 Translation Error
Error occurred while flattening model HelmholtzMedia.Examples.BranchingDynamicPipes

reference state is fixed

Currently, the reference state is fixed.
To change the reference state, one has to change the helmholtzCoefficients.idealPower.
Examples.Validation.ReferenceState contains some more explanation

convergence criterion

Revise convergence criteria for all iterative procedures.
Consider not just residual, but also size of Newton step.
Especially continue iterating as long as density still sees updates.

go parallel

I don't see how to do this in Modelica, but going parallel could significantly speed up things:
All 5-14 HelmholtzDerivs can be calculated simultaneously, see HelmholtzDerivs and setHelmholtzDerivs
and all 12 up to 50 terms of each HelmholtzDeriv could be evaluated simultaneously, see e.g. f_r

Combining these two should give 60 or more independent threads.
This should be investigated in combination with (automatic?) common subexpression elimination, because the f_r etc function calls have very many terms in common! Maybe combine with #26 ?

Helium: convergence of setState_ph and setState_ps

For Helium, the two functions setState_ph and setState_ph have converge problems in the super-critical liquid-like region and in the super-critical vapour-like region.
Three possible solutions:

• Better start values. For water there the IAPWS backward equations, similar equations could be set up for other fluids. But that means quite a lot of work. I could not find any other, more general approaches for guessing start values.
• Use a iterative procedure other than Newton. For example, Broyden (2D Secant) or nested 1D Ridders/Bisection/Brent/Regula falsi.
• Make Newton more stable using linesearching and backtracking,
  as described by e.g. Numerical Recipes 3rd ed.
  http://books.google.de/books?id=DyykEZo4fwUC&pg=PA477
  or by Numerik für Ingenieure und Naturwissenschaftler, 2nd ed., Dahmen (2008)
  http://books.google.de/books?id=zWenT-hxDxEC&pg=PA200

Problems with Helium model

We are trying to use the Helium model for cryogenic system studies, possibly using OpenModelica. We already found some issues using Dymola. Please download this small package which allows you to reproduce them: http://home.deib.polimi.it/casella/transfer/HelmholtzBugs.mo

We have built a test case, named CheckDerivativesHeliumHelmholzFluids, to check the consistency of density derivatives by integrating them, and then comparing the result with the actual density changes. The model has been tested with other media (e.g., IF97 water, or ExternalMedia models of Helium) and it works. When using the Helium model of your library, if fails with this error message:

The following error was detected at time: 3.901681204235383E-005
Out of memory for array dimensions
It could due to too many matrices, infinite recursion, or uninitialized variables.
You can increase the size of 'Sizebuffer' in dymola/source/matrixop.h.
The stack of functions is:
HelmholtzMedia.HelmholtzFluids.Helium.EoS.f_i_Unique9
HelmholtzMedia.HelmholtzFluids.Helium.EoS.setHelmholtzDerivsSecond_Unique8
HelmholtzMedia.HelmholtzFluids.Helium.setState_phX_Unique42
HelmholtzMedia.Interfaces.PartialHelmholtzMedium.setState_ph_Unique41
HelmholtzMedia.Interfaces.PartialHelmholtzMedium.setState_ph_Unique41(p0, h, 0)

It looks pretty much low-level, can you have a look at it?