Details on PropertyTester

[desyk]

Hello all,

I have developed a COBIA Standalone Thermodynamic Property Package which uses a 3rd-order Central Difference Method for the calculation of the various derivatives, as described in this article. It works quite well, producing results that are within 0.1% of those calculated trough perturbation by the PropertyTester Unit Operation (with few exceptions when the derivative is very close to zero).
Out of curiosity, and as I believe the latter are more accurate, could anyone share details on the numerical method PropertyTester uses to calculate the derivatives?

Best,
D.S.

[jasper]

Temperature and pressure derivatives are perturbed by a relative distance that can be adjusted from the unit operation's configuration; the perturbation is in the direction of positive T and P.

Mole fraction derivatives are unconstrained and perturbed with a perturbation relative to 1, so the relative perturbation is the absolute perturbation for mole fractions. For X > 0.5, the direction is negative, otherwise positive. These derivatives should be taken with a grain of salt, as it is possible that any component may re-normalize the compositions, throwing off the perturbations.

Mole number derivatives are similar to mole fraction derivatives, except the perturbed composition is normalized. For extensive properties, e.g. molar enthalpy h, the derivative is for the total enthalpy H = n * h, where n is the total mole number, and the derivative is defined for a total number of moles of n = 1. So in this case chain rule must be applied and the derivative must be corrected for h itself. So effectively the perturbed property enthalpy.Dmoles is h(T,P,x) + (h(T,P,x1) - h(T,P,x))/delta where x1 is offset by delta for the component that is perturned, and then normalized: x1 = (x + [0,0,...,delta,...,0])/(1+delta)

[desyk]

Thank you for the information. I went ahead and implemented the mole number derivatives (I had only implemented Temperature and Pressure derivatives so far). The only difference to the method you described is that each derivative is calculated using the Central Finite Difference method. This works very well for the Mole number derivatives of intensive properties, but fails when it comes to extensive properties.
I suspect this is because I am using REFPROP to retrieve those properties directly in their extensive form -- i.e. I do not calculate the total enthalpy in order to get its derivative. Instead, I perform the last step you described from the get-go, as such:

enthalpy.Dmoles = [(1/12) * h(T,P,x1) + (-2/3) * h(T,P,x2) + (2/3) * h(T,P,x3) + (-1/12) * h(T,P,x4)] / dx, where:
x1 = x0 - [0,0, ..., 2*dx, ..., 0] / (1 - 2*dx)
x2 = x0 - [0,0, ..., dx, ..., 0] / (1 - dx)
x3 = x0 + [0,0, ..., dx, ..., 0] / (1 + dx)
x4 = x0 + [0,0, ..., 2*dx, ..., 0] / (1 + 2*dx)

This calculation yields:
[12444.8, -2207.25, -2144.26]
Vs. the perturbed:
[26289.1, 11638.1, 11701.6]
(note the direction change)

I also performed some simple manual testing by inputting the x0..4 composition vectors straight into REFPROP and calculating the derivatives in a spreadsheet using your method, but got fairly similar results to the ones given by my implementation.

The same mixture in TEA yields:
[3965.16, 2969.46, 3022.72] (TEA)
VS
[3965.08, 2969.44, 3022.73] (Perturbed)
What I find interesting (and confusing) is that it reports a positive derivative for both the 2nd and 3rd compounds, but, when I manually alter the composition so that i.e. the molar fraction of the 2nd compound increases, the molar entalpy of the mixture decreases. Considering that, a negative derivative would make sense in this case. Did I miss something from the Thermo Specification which would explain this behaviour?


The molar fractions of the mixture tested were:
x[CO2] = 0.15
x[N2] = 0.7
x[O2] = 0.15