github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-04-01 03:00:13 -04:00
Code Issues Packages Projects Releases Wiki Activity

Added code for liquid enthalpy ancillary for Propane - proof of principle.

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-07-31 19:39:24 +02:00
parent 33b5fe0e28
commit fc7eb55ca2
8 changed files with 174 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
dev/fluids/n-Propane.json
View File

@@ -127,6 +127,11 @@
"type": "rhoV",
"using_tau_r": true
},
"hL":{
"A": [-6.0513675952483362e-13, 8.8373829521070838e-10, -5.7556127200695971e-07, 0.00020526029744824865, -0.042961977764103827, 0.5027557264334126, 2471.9788731127519, -348676.54635279981],
"B":[-0.0026659218650113606, 1],
"type" : "rational_polynomial"
},
"surface_tension": {
"BibTeX": "Mulero-JPCRD-2012",
"Tc": 369.89,

18
dev/scripts/fit_rational_functions.py
View File

@@ -42,26 +42,28 @@ def fit_rational_polynomial(x, y, xfine, n, d):
AB = scipy.optimize.curve_fit(obj, x, y, p0 = ABlin)[0]
poles = np.roots(list(ABnonlin[n+1::])+[1])
poles = np.roots(list(AB[n+1::])+[1])
poles = poles[np.isreal(poles)]
poles = poles[poles > min(x)]
poles = poles[poles < max(x)]
else:
# Find a pole that is not in the range of x
def obj2(Tpole,x,y,AB):
def obj2(Tpole, x, y, AB):
B = -1/Tpole
A = np.polyfit(x,y*(x*B+1),n)
yfit = np.polyval(A,x)/(x*B + 1)
A = np.polyfit(x, y*(x*B+1), n)
yfit = np.polyval(A, x)/(x*B + 1)
AB[:] = list(A) + [B] # Set so that it uses the AB passed in rather than making local variable
rms = np.sqrt(np.sum(np.power(yfit-y,2)))
rms = np.sqrt(np.sum(np.power(yfit-y, 2)))
return rms
AB = []
scipy.optimize.fminbound(obj2, Tc+0.1, 1.5*Tc, args = (x, y, AB))
return dict(max_abs_error = np.max(np.abs(obj(x, *AB)-y)),
yfitnonlin = obj(xfine, *AB)
yfitnonlin = obj(xfine, *AB),
A = AB[0:n+1],
B = list(AB[n+1::]) + [1]
)
class SplineFitter(object):
@@ -115,7 +117,7 @@ def strictly_decreasing(L):
from matplotlib.backends.backend_pdf import PdfPages
pp = PdfPages('multipage.pdf')
for i, fluid in enumerate(sorted(CoolProp.__fluids__)):
for i, fluid in enumerate(['n-Propane']):#sorted(CoolProp.__fluids__)):
fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2, 2)
@@ -151,6 +153,8 @@ for i, fluid in enumerate(sorted(CoolProp.__fluids__)):
ax1.plot(xfine, rp['yfitnonlin'] + rp['max_abs_error'], 'k--')
ax1.plot(xfine, rp['yfitnonlin'] - rp['max_abs_error'], 'k--')
print rp['A'], rp['B']
rp = fit_rational_polynomial(x, hfg, xfine, n, d)
ax2.plot(x, hfg)
ax2.plot(xfine, rp['yfitnonlin'],'r')

88
include/CoolPropFluid.h
View File

@@ -18,6 +18,8 @@
#include <map>
#include <assert.h>
#include <iterator>
#include "Eigen/Core"
#include "PolyMath.h"
namespace CoolProp {
@@ -289,55 +291,89 @@ public:
class SaturationAncillaryFunction
{
private:
Eigen::MatrixXd num_coeffs, ///< Coefficients for numerator in rational polynomial
den_coeffs; ///< Coefficients for denominator in rational polynomial
std::vector<double> n, t, s;
bool using_tau_r;
double Tmax, Tmin, reducing_value, T_r;
long double Tmax, Tmin, reducing_value, T_r, max_abs_error;
int type;
enum ancillaryfunctiontypes{TYPE_NOT_EXPONENTIAL = 0, TYPE_EXPONENTIAL = 1};
enum ancillaryfunctiontypes{TYPE_NOT_SET = -1, TYPE_NOT_EXPONENTIAL = 0, TYPE_EXPONENTIAL = 1, TYPE_RATIONAL_POLYNOMIAL = 2};
std::size_t N;
public:
SaturationAncillaryFunction(){};
SaturationAncillaryFunction(){type = TYPE_NOT_SET;};
SaturationAncillaryFunction(rapidjson::Value &json_code)
{
n = cpjson::get_double_array(json_code["n"]);
t = cpjson::get_double_array(json_code["t"]);
Tmin = cpjson::get_double(json_code,"Tmin");
Tmax = cpjson::get_double(json_code,"Tmax");
reducing_value = cpjson::get_double(json_code,"reducing_value");
using_tau_r = cpjson::get_bool(json_code,"using_tau_r");
T_r = cpjson::get_double(json_code,"T_r");
std::string type = cpjson::get_string(json_code,"type");
if (!type.compare("rhoLnoexp"))
if (!type.compare("rational_polynomial"))
{
std::vector<double> A = cpjson::get_double_array(json_code["A"]);
std::vector<double> B = cpjson::get_double_array(json_code["B"]);
std::reverse(A.begin(), A.end());
std::reverse(B.begin(), B.end());
num_coeffs = vec_to_eigen(A);
den_coeffs = vec_to_eigen(B);
}
else
{
n = cpjson::get_double_array(json_code["n"]);
t = cpjson::get_double_array(json_code["t"]);
Tmin = cpjson::get_double(json_code,"Tmin");
Tmax = cpjson::get_double(json_code,"Tmax");
reducing_value = cpjson::get_double(json_code,"reducing_value");
using_tau_r = cpjson::get_bool(json_code,"using_tau_r");
T_r = cpjson::get_double(json_code,"T_r");
}
if (!type.compare("rational_polynomial"))
this->type = TYPE_RATIONAL_POLYNOMIAL;
else if (!type.compare("rhoLnoexp"))
this->type = TYPE_NOT_EXPONENTIAL;
else
this->type = TYPE_EXPONENTIAL;
this->N = n.size();
s = n;
};
/// Get the maximum absolute error for this fit
long double get_max_abs_error(){return max_abs_error;};
double evaluate(double T)
{
double THETA = 1-T/T_r;
for (std::size_t i = 0; i < N; ++i)
if (type == TYPE_NOT_SET)
{
s[i] = n[i]*pow(THETA, t[i]);
throw ValueError(format("type not set"));
}
double summer = std::accumulate(s.begin(), s.end(), 0.0);
if (type == TYPE_NOT_EXPONENTIAL)
else if (type == TYPE_RATIONAL_POLYNOMIAL)
{
return reducing_value*(1+summer);
Polynomial2D poly;
return poly.evaluate(num_coeffs, T)/poly.evaluate(den_coeffs, T);
}
else
{
double tau_r_value;
if (using_tau_r)
tau_r_value = T_r/T;
double THETA = 1-T/T_r;
for (std::size_t i = 0; i < N; ++i)
{
s[i] = n[i]*pow(THETA, t[i]);
}
double summer = std::accumulate(s.begin(), s.end(), 0.0);
if (type == TYPE_NOT_EXPONENTIAL)
{
return reducing_value*(1+summer);
}
else
tau_r_value = 1.0;
return reducing_value*exp(tau_r_value*summer);
{
double tau_r_value;
if (using_tau_r)
tau_r_value = T_r/T;
else
tau_r_value = 1.0;
return reducing_value*exp(tau_r_value*summer);
}
}
}
double invert(double value)
@@ -464,7 +500,7 @@ public:
struct Ancillaries
{
SaturationAncillaryFunction pL, pV, rhoL, rhoV;
SaturationAncillaryFunction pL, pV, rhoL, rhoV, hL, hV, sL, sV;
MeltingLineVariables melting_line;
SurfaceTensionCorrelation surface_tension;
};

6
include/PolyMath.h
View File

@@ -258,7 +258,7 @@ public:
/// @param coefficients vector containing the ordered coefficients
/// @param x_in double value that represents the current input in the 1st dimension
/// @param firstExponent integer value that represents the lowest exponent of the polynomial
/// @param x_base double value that represents the base value for a centred fit in the 1st dimension
/// @param x_base double value that represents the base value for a centered fit in the 1st dimension
double evaluate(const Eigen::MatrixXd &coefficients, const double &x_in, const int &firstExponent=0, const double &x_base=0.0);
/// @param coefficients matrix containing the ordered coefficients
@@ -266,8 +266,8 @@ public:
/// @param y_in double value that represents the current input in the 2nd dimension
/// @param x_exp integer value that represents the lowest exponent of the polynomial in the 1st dimension
/// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension
/// @param x_base double value that represents the base value for a centred fit in the 1st dimension
/// @param y_base double value that represents the base value for a centred fit in the 2nd dimension
/// @param x_base double value that represents the base value for a centered fit in the 1st dimension
/// @param y_base double value that represents the base value for a centered fit in the 2nd dimension
double evaluate(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &x_exp, const int &y_exp, const double &x_base=0.0, const double &y_base=0.0);
/// @param coefficients vector containing the ordered coefficients

54
src/Backends/Helmholtz/FlashRoutines.cpp
View File

@@ -53,7 +53,7 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
rhoLsat = HEOS.solver_rho_Tp(HEOS._T, psatLanc, rhoLanc);
rhoVsat = HEOS.solver_rho_Tp(HEOS._T, psatVanc, rhoVanc);
if (!ValidNumber(rhoLsat) || !ValidNumber(rhoVsat) ||
fabs(rhoLsat/rhoLanc-1) > 0.1 || fabs(rhoVanc/rhoVsat-1) > 0.1)
fabs(rhoLsat/rhoLanc-1) > 0.5 || fabs(rhoVanc/rhoVsat-1) > 0.5)
{
throw ValueError("pseudo-pure failed");
}
@@ -325,9 +325,59 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
throw NotImplementedError("PHSU_D_flash not ready for mixtures");
}
}
void FlashRoutines::HSU_P_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, int other, long double T0, long double rhomolar0)
{
}
void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
{
throw NotImplementedError("HSU_P_flash Not implemented yet");
if (HEOS.imposed_phase_index > -1)
{
// Use the phase defined by the imposed phase
HEOS._phase = HEOS.imposed_phase_index;
}
else
{
if (HEOS.is_pure_or_pseudopure)
{
// Find the phase, while updating all internal variables possible
switch (other)
{
case iSmolar:
HEOS.p_phase_determination_pure_or_pseudopure(iSmolar, HEOS._smolar); break;
case iHmolar:
HEOS.p_phase_determination_pure_or_pseudopure(iHmolar, HEOS._hmolar); break;
case iUmolar:
HEOS.p_phase_determination_pure_or_pseudopure(iUmolar, HEOS._umolar); break;
default:
throw ValueError(format("Input for other [%s] is invalid", get_parameter_information(other, "long").c_str()));
}
}
else
{
HEOS._phase = iphase_gas;
throw NotImplementedError("HSU_P_flash does not support mixtures (yet)");
// Find the phase, while updating all internal variables possible
}
}
if (HEOS.isHomogeneousPhase() && !ValidNumber(HEOS._p))
{
switch (other)
{
case iDmolar:
break;
case iHmolar:
HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._hmolar, iHmolar); break;
case iSmolar:
HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._smolar, iSmolar); break;
case iUmolar:
HEOS._rhomolar = HEOS.solver_for_rho_given_T_oneof_HSU(HEOS._T, HEOS._umolar, iUmolar); break;
default:
break;
}
HEOS.calc_pressure();
HEOS._Q = -1;
}
}
void FlashRoutines::DHSU_T_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
{

9
src/Backends/Helmholtz/FlashRoutines.h
View File

@@ -49,6 +49,13 @@ public:
/// @param other The index for the other input from CoolProp::parameters; allowed values are iHmolar, iSmolar, iUmolar
static void HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other);
/// The single-phase flash routine for the pairs (P,H), (P,S), and (P,U). Similar analysis is needed
/// @param HEOS The HelmholtzEOSMixtureBackend to be used
/// @param other The index for the other input from CoolProp::parameters; allowed values are iHmolar, iSmolar, iUmolar
/// @param T0 The initial guess value for the temperature [K]
/// @param rho0 The initial guess value for the density [mol/m^3]
static void HSU_P_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, int other, long double T0, long double rhomolar0);
/// A generic flash routine for the pairs (D,P), (D,H), (D,S), and (D,U). Similar analysis is needed
/// @param HEOS The HelmholtzEOSMixtureBackend to be used
/// @param other The index for the other input from CoolProp::parameters; allowed values are iP, iHmolar, iSmolar, iUmolar
@@ -56,4 +63,4 @@ public:
};
} /* namespace CoolProp */
#endif /* FLASHROUTINES_H */
#endif /* FLASHROUTINES_H */

10
src/Backends/Helmholtz/Fluids/FluidLibrary.h
View File

@@ -926,6 +926,16 @@ protected:
fluid.ancillaries.pV = SaturationAncillaryFunction(ancillaries["pV"]);
fluid.ancillaries.rhoL = SaturationAncillaryFunction(ancillaries["rhoL"]);
fluid.ancillaries.rhoV = SaturationAncillaryFunction(ancillaries["rhoV"]);
if (ancillaries.HasMember("hL"))
{
fluid.ancillaries.hL = SaturationAncillaryFunction(ancillaries["hL"]);
}
else
{
if (get_debug_level() > 0){ std::cout << "Missing hL ancillary for fluid " << fluid.name; }
}
};
/// Parse the surface_tension

23
src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
View File

@@ -625,6 +625,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot
switch (other)
{
// First try the ancillaries, use them to determine the state if you can
case iT:
{
long double p_vap = 0.98*static_cast<double>(_pVanc);
@@ -638,6 +639,28 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot
}
break;
}
case iHmolar:
{
long double h_liq = components[0]->ancillaries.hL.evaluate(_TLanc);
long double h_vap = components[0]->ancillaries.hV.evaluate(_T);
// Check if in range given the accuracy of the fit
if (value > h_vap + components[0]->ancillaries.hV.get_max_abs_error()){
this->_phase = iphase_gas; _Q = -1000; return;
}
else if (value < h_liq - components[0]->ancillaries.hL.get_max_abs_error()){
this->_phase = iphase_liquid; _Q = 1000; return;
}
break;
}
case iSmolar:
{
// Add entropy ancillary code here
}
case iUmolar:
{
// Add entropy ancillary code here
}
default:
{
// Always calculate the densities using the ancillaries