Implementing tabular heat capacity
 

Hello,

I'm currently working on implementing a heat capacity model using tabulated T- and p-dependent data.

I based the library on the hPolynomial model which works with my solver, and reimplemented cp(), ha(), hs(), hc() and s() using an interpolation2DTable for cp and a T-integrated table for enthalpy calculation.

The values in the table are comparable to the values I used in the polynomial model, however, the temperature calculation fails on the very first step, with the "Maximum number of iterations exceeded" error in the temperature calculation function T().

When I use the limit() function to specify valid ranges for T, I see very extreme negative values (-90000 to -1e7). Clamping them does not help. Why would there be negative temperature values? Are they an artifact of the iteration procedure in T()? Is it diverging because I can't provide meaningful data in this range? I have the tables set to clamp, so the capacity and the enthalpy will be constant for these temperatures. Would it help if I use linear extrapolation instead of clamping maybe?

Is there anything tricky that I might have missed? The tabulated data isn't as smooth because it contains enthalpies of transitions, can this lead to problems?

I'm also a bit confused about the mixture functionality. I only want to work with a single component, does this mean I can safely ignore the operatores for adding, subtracting and multiplying mixtures? I don't have meaningful implementations of them right now, but to my current understanding it should not be needed for my case?

On a related note, where can I find information about the temperature calculation? I'm having some difficulties understanding the code in the T() function and it's subroutines.
I'm also not quite sure where the reference temperature specie::Tstd comes from? I have looked into constants.h and its includes, but where is the value specified?

Edit: I have looked some more into the temperature iteration routine, and I believe it uses a newton algorithm for inverting the H(T,p) function to T. This points indeed to problems when the gradient of this function becomes zero, as it does outside of the table range. I'll try to implement an extrapolating table class to fix this. However, please share your knowledge with me ;)

I have created an extrapolating2DTable class which I successfully tested with some simple values, however I'm still getting this error in the temperature calculation.

I guess you tried this already but what happens if your heat capacity values are a constant over the full temperature range?

So far I have not yet had much success, although some very simple cases seem to work (such as constant heat capacity and constant energy in the case so the temperature doesn't change).

I have tried to decrease the step in the temperature approximation newton method with some limited success, at the cost of having to use more steps until it converges. Using that, I was able to use a linear increasing heat capacity. I think that maybe the data I'm using was too rough, I'll have to check that.

For example: Linear increasing heat capacity using two data points. Enthalpy should then be quadratic in temperature, however you can't have that with only two points.

I'll try to see if this is one of the problems, however, I was also using finer data which also was problematic.

I believe that something else is flawed, as the Newton method shouldn't need so many steps with linear data.

The heat capacity is strictly positive, so there should not be a problem with Newtons method AFAIK.

I finally got it working, turns out I had incorrect enthalpy tables.

I stumbled upon further problems, which are very likely caused by the tables not being continously differentiable, which is a prerequisite for the Newton method for root-finding.

Now I'm wondering, is there a more robust method for finding roots that doesn't require continous differentiable functions, or do I need to implement an extrapolating table that employs spline interpolation?

Edit: I implemented a bisection method in addition to the newton method used for inverting f(T). On the first look it appears to work, but I'm having other (mostly unrelated) problems which might interfere, so I need to do further comparisons.