search

Cantera / cantera

Chemical kinetics, thermodynamics, and transport tool suite
https://cantera.org
Other
604 stars 346 forks source link

Incorrect heat capacities near liquid/vapor transition for "TPX" species #273

Closed speth closed 7 years ago

speth commented 9 years ago

Calculating specific heat capacities (cp or cv) for states near the saturation curve give incorrect results. The values appear to asymptotically approach the values which are returned inside the liquid/vapor region (where the reported values of cp are incorrect -- cp should be infinite there). For example:

import cantera as ct
g = ct.Nitrogen()
for P in [760000, 770000, 780000, 780100, 780400, 780450, 780455, 
          780500, 780600, 781000, 782000, 790000]:
    g.TP = 100, P
    print P, g.X, g.cv, g.cp

which gives:

760000 1.0 838.402563772 1454.45597099
770000 1.0 840.199857557 1464.58564952
780000 1.0 842.018464973 8029648.73494
780100 1.0 1592.01729477 8029371.99696
780400 1.0 4848.92500893 8028541.90605
780450 1.0 5391.40302329 8028403.54613
780455 1.0 5445.69230038 8028389.70217
780500 0.0 1547.70533543 8028265.19516
780600 0.0 1541.17494617 8027988.53878
781000 0.0 1515.0538804 2310.60053871
782000 0.0 1449.75160413 2310.53856077
790000 0.0 935.06668906 2310.05118594

In comparison to the results which can be obtained from the NIST Chemistry Webbook.

This behavior is consistent across all of the "TPX" species, so this is an issue with its algorithm for evaluating heat capacities, rather than an issue with any of the particular equations of state.

bryanwweber commented 8 years ago

This is because the specific heats are evaluated numerically as T*(s(T+dt) - s(T-dt))/((T+dt) - (T-dt)) at either constant volume or constant pressure. If T is close enough to T_sat, adding or subtracting dT will push one point or the other into the vapor dome, resulting in the incorrect values shown here. It seems a new algorithm for calculating the derivative is necessary.

speth commented 8 years ago

Calculating cv and cp by finite difference seems a bit silly when s is computed by integrating the expression for cv(T) in the first place. But implementing that would require a separate method for each species, rather than a single implementation in Substance, which is probably why it was initially done this way.

bryanwweber commented 8 years ago

Only part of s for the TPX substances is calculated by integrating cv(T). In particular, the only portion of cv(T) given by the curve fits in Reynolds' book is the zero-density limit, so calculating cv(T,p) for substances not near that limit is not possible with only the given coefficients (the G array). I don't think writing an implementation in each Substance is a particular disadvantage, provided that each substance had an accurate, explicit formula - the trouble is, they don't (as far as I know). I suppose it would hypothetically be possible to analytically differentiate the formula for s for each Substance, but that sounds like a big pain, given the complexity of the EOS.

It seems like it should be fairly easy to improve this situation by either ensuring that the step size in T for the current formulation doesn't go into the two-phase region, using a higher-order differentiation, possibly with varying step sizes, or both. This would likely come at the cost of increased computational time though.