ThermoState 0.4.5: specification of thermodynamic properties

Hi everybody. i’m been working on this package for already some time, and it seems like is the time to make an announcement on discourse

What does it do?

if i oversymplify, is an equivalent of Tables.jl for Thermodynamic Models.

A problem when designing those models is that you want a lot of properties with different implementations, For example, if you want the the enthalpy of a liquid, maybe you want to calculate that using Pressure and Temperature, or Pressure and Entropy, or Molar Volume and Temperature, or specifify that is saturated… etc. at the start, if you have only one model, this is easy, you can just write the implementations directly:

pressure(v,t)
entropy(v,t)
...

if you have two or more models, maybe do something like this:

pressure(model::Model1,v,t)
entropy(model::Model1,v,t)
...
pressure(model::Model2,v,t)
entropy(model::Model2,v,t)

the problem arrises when now you want to add a different way to calculate said properties. one way is to add a function preffix or suffix:

#volume - temperature
pressure_vt(model::Model1,v,t)
entropy_vt(model::Model1,v,t)
...
pressure_vt(model::Model2,v,t)
entropy_vt(model::Model2,v,t)

#entropy - temperature
pressure_st(model::Model1,s,t)
entropy_st(model::Model1,s,t)
...
pressure_vt(model::Model2,s,t)
entropy_vt(model::Model2,s,t)

other option is to use a type system to signal on what function dispatch

#volume - temperature
pressure_vt(model::Model1,v,t)
entropy_vt(model::Model1,v,t)
...
pressure(model::Model2,::VT,v,t)
entropy(model::Model2,::VT,v,t)

#entropy - temperature
pressure(model::Model1,::ST,s,t)
entropy(model::Model1,::ST,s,t)
...
pressure_vt(model::Model2,::ST,s,t)
entropy_vt(model::Model2,::ST,s,t)

either way, the user interface has become cluttered. with multicomponent models, the problem is aggravated, as now the material parts can come in diferent forms, do i want a mol fraction, a mass fraction? total amounts? .

With all done, now there exists the Unitful.jl package, that, while excelent on what it does, now requires another dispatch.

This package was done to explore a possible solution of this problem.

Usage:

The package revolves around the ThermodynamicState object, that can be created with the state function:

st = state(t=1.0u"°C",p=1.0u"atm")

all property functions have the form : property(model,state,units,args..,kwargs...) the “Simplest model” is the one that doesn’t calculate anything, and returns the values stored in thermodynamic state:

t0 = temperature(FromState(),st) #kelvin by default, no Unitful, returns 274.15
t0 = temperature(FromState(),st,u"°C") #no Unitful, returns 1.0
t0 = @to_units temperature(FromState(),st) #kelvin by default, with Unitful, returns 274.15 K
t0 = @to_units temperature(FromState(),st,u"°C") #with Unitful, returns 1.0 °C

there is a lot more to explain, in the Readme is a somewhat more detailed explanation of what it does.

I have been eating my own dog food and using this package to port some thermodynamic models to the interface proposed here. some resulting packages are:
(almost registered, properties of moist air)


(not registered, collection of models)

If this is useful for someone else besides myself, let me know if i can help.

4 Likes

Why did you choose to return unitless values by default then add a macro to get units? Seems backwards to me. There are already functions in unit full to strip off units.

Seems like a good goal and approach, I’ve run into lesser versions of the same problem

1 Like

At the start, i did return Unitful Units, but i run of some problems. for example. in a Antoine implementation:
p = exp(a + b/(c+t))
this would require that the implementation itself strips the numbers,or add units to every constant to make the expression inside the exp adimensional. passing just numbers allows passing the properties to the implementations easier. Also, an Unitful Quantity is a number, not a Real, and had blocked me from using ForwardDiff. by passing numbers, an hypothetical user can do this:

∂p∂t(t) = ForwardDiff.derivative(tt->pressure(model,st(tt)),t) #using variable specifications

this is somewhat alleviated in the fact that when printing an state that is not normalized, the units are shown:

julia> st = state(t=1u"°C",p=2u"atm",normalize_units=false)
ThermodynamicState with 2 properties:
  Temperature : 1[°C]
  Pressure : 2[atm]

But isn’t this the idea?

2 Likes

If you are manipulating the variables in a physically consistent system, the units help a lot . but a lot of property models are … blackboxes, :smile:. , non consistent equations, magical constants. Even the authors of said models didnt add any units to their models, as those models are the result of an optimization algorithm.
Returning units from the property functions can be be solved with adding one uconvert(preferred.unit,val) to each value passed to each model in each implementation, and a second uconvert from the model’s output models to the desired output units. so its doable. it just transfers the responsability of converting from the state to who implements the model. i decided that the state has to have the responsability of managing the unit conversions.

I was talking about this pretty extensively with Tim Holy the other day: Plotting Unitful Functions

We basically have to fix all of science to do this properly :sweat_smile:

4 Likes

You are manipulating units in a physically consistent system, otherwise it wouldn’t work. Its just that some people chose to leave off their units, which is basically what Raf’s link says I think. All the constants in those equations have units, which you can derive based on statements like “this argument should be in kelvin”.

So an alternate path would be to derive the units for the constants, put them in, then define something like
exp(x::Unitful.Quantity) = exp(uconvert(Unitful.NoUnits, x))

I think the hardest cases are things like P(T, n) = A*T^n because the units of where T is a temperature and n is a unitless exponent that may be a fit parameter. In that case the units of A depend on n. edit: In the other thread they suggested doing P(T, n) = A*(P/1u"Torr")^n for cases like this, which allows A to have fixed units.

1 Like