github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-02-09 05:15:45 -05:00
Code Issues Packages Projects Releases Wiki Activity

Lots more work on melting curves, can now make pretty phase diagram for water including the melting line

Browse Source
This commit is contained in:
Ian Bell
2014-08-10 19:28:04 +02:00
parent 51b34dcdbe
commit 36d601490f
12 changed files with 160 additions and 50 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
src/AbstractState.cpp
View File

@@ -186,6 +186,11 @@ double AbstractState::trivial_keyed_output(int key)
double AbstractState::keyed_output(int key)
{
if (get_debug_level()>=50) std::cout << format("AbstractState: keyed_output called for %s ",get_parameter_information(key,"short").c_str()) << std::endl;
// Handle trivial inputs
if (is_trivial_parameter(key))
{
return trivial_keyed_output(key);
}
switch (key)
{
case iQ:

123
src/Backends/Helmholtz/Fluids/Ancillaries.cpp
View File

@@ -1,5 +1,13 @@
#include "Ancillaries.h"
#include "DataStructures.h"
#include "AbstractState.h"
#if defined(ENABLE_CATCH)
#include "crossplatform_shared_ptr.h"
#include "catch.hpp"
#endif
namespace CoolProp{
@@ -103,20 +111,29 @@ double SaturationAncillaryFunction::invert(double value)
void MeltingLineVariables::set_limits(void)
{
if (type == MELTING_LINE_SIMON_TYPE){
MeltingLinePiecewiseSimonSegment &partmin = simon.parts[0];
MeltingLinePiecewiseSimonSegment &partmax = simon.parts[simon.parts.size()-1];
Tmin = partmin.T_0;
Tmax = partmax.T_max;
pmin = partmin.p_0;
pmax = evaluate(iP, iT, Tmax);
// Fill in the min and max pressures for each part
for (std::size_t i = 0; i < simon.parts.size(); ++i){
MeltingLinePiecewiseSimonSegment &part = simon.parts[i];
part.p_min = part.p_0 + part.a*(pow(part.T_min/part.T_0,part.c)-1);
part.p_max = part.p_0 + part.a*(pow(part.T_max/part.T_0,part.c)-1);
}
pmin = simon.parts.front().p_min;
pmax = simon.parts.back().p_max;
Tmin = simon.parts.front().T_min;
Tmax = simon.parts.back().T_max;
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE){
MeltingLinePiecewisePolynomialInTrSegment &partmin = polynomial_in_Tr.parts[0];
MeltingLinePiecewisePolynomialInTrSegment &partmax = polynomial_in_Tr.parts[polynomial_in_Tr.parts.size() - 1];
Tmin = partmin.T_0;
Tmax = partmax.T_max;
pmin = partmin.p_0;
pmax = evaluate(iP, iT, Tmax);
// Fill in the min and max pressures for each part
for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i){
MeltingLinePiecewisePolynomialInTrSegment &part = polynomial_in_Tr.parts[i];
part.p_min = part.evaluate(part.T_min);
part.p_max = part.evaluate(part.T_max);
}
Tmin = polynomial_in_Tr.parts.front().T_min;
pmin = polynomial_in_Tr.parts.front().p_min;
Tmax = polynomial_in_Tr.parts.back().T_max;
pmax = polynomial_in_Tr.parts.back().p_max;
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE){
MeltingLinePiecewisePolynomialInThetaSegment &partmin = polynomial_in_Theta.parts[0];
@@ -124,7 +141,7 @@ void MeltingLineVariables::set_limits(void)
Tmin = partmin.T_0;
Tmax = partmax.T_max;
pmin = partmin.p_0;
pmax = evaluate(iP, iT, Tmax);
//pmax = evaluate(iP, iT, Tmax);
}
else{
throw ValueError("only Simon supported now");
@@ -151,11 +168,7 @@ long double MeltingLineVariables::evaluate(int OF, int GIVEN, long double value)
for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i){
MeltingLinePiecewisePolynomialInTrSegment &part = polynomial_in_Tr.parts[i];
if (is_in_closed_range(part.T_min, part.T_max, T)){
long double summer = 0;
for (std::size_t i =0; i < part.a.size(); ++i){
summer += part.a[i]*(pow(T/part.T_0,part.t[i])-1);
}
return part.p_0*(1+summer);
return part.evaluate(T);
}
}
throw ValueError("unable to calculate melting line (p,T) for polynomial_in_Tr curve");
@@ -189,23 +202,21 @@ long double MeltingLineVariables::evaluate(int OF, int GIVEN, long double value)
return T;
}
}
throw ValueError("unable to calculate melting line (p,T) for Simon curve");
throw ValueError("unable to calculate melting line p(T) for Simon curve");
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE || type == MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE)
else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE)
{
class solver_resid : public FuncWrapper1D
class solver_resid : public FuncWrapper1D
{
public:
MeltingLineVariables *line;
MeltingLinePiecewisePolynomialInTrSegment *part;
long double r, given_p, calc_p, T;
solver_resid(MeltingLineVariables *line, long double p) : line(line), given_p(p){};
solver_resid(MeltingLinePiecewisePolynomialInTrSegment *part, long double p) : part(part), given_p(p){};
double call(double T){
this->T = T;
// Calculate p using melting line
calc_p = line->evaluate(iP, iT, T);
calc_p = part->evaluate(T);
// Difference between the two is to be driven to zero
r = given_p - calc_p;
@@ -213,17 +224,63 @@ long double MeltingLineVariables::evaluate(int OF, int GIVEN, long double value)
return r;
};
};
solver_resid resid(this, value);
double pmin = evaluate(iP, iT, Tmin);
double pmax = evaluate(iP, iT, Tmax);
std::string errstr;
return Brent(resid, Tmin, Tmax, DBL_EPSILON, 1e-12, 100, errstr);
// Need to find the right segment
for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i){
MeltingLinePiecewisePolynomialInTrSegment &part = polynomial_in_Tr.parts[i];
if (is_in_closed_range(part.p_min, part.p_max, value)){
std::string errstr;
solver_resid resid(&part, value);
double T = Brent(resid, part.T_min, part.T_max, DBL_EPSILON, 1e-12, 100, errstr);
return T;
}
}
throw ValueError("unable to calculate melting line T(p) for polynomial_in_Tr curve");
}
else{
throw ValueError(format("Invalid melting line type (T,p) [%d]",type));
throw ValueError(format("Invalid melting line type T(p) [%d]",type));
}
}
}
}; /* namespace CoolProp */
#if defined(ENABLE_CATCH)
TEST_CASE("Water melting line", "")
{
shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS","water"));
int iT = CoolProp::iT, iP = CoolProp::iP;
SECTION("Ice Ih-liquid")
{
double actual = AS->melting_line(iT, iP, 138.268e6);
double expected = 260.0;
CAPTURE(actual);
CAPTURE(expected);
CHECK(std::abs(actual-expected) < 0.01);
}
SECTION("Ice III-liquid")
{
double actual = AS->melting_line(iT, iP, 268.685e6);
double expected = 254;
CAPTURE(actual);
CAPTURE(expected);
CHECK(std::abs(actual-expected) < 0.01);
}
SECTION("Ice V-liquid")
{
double actual = AS->melting_line(iT, iP, 479.640e6);
double expected = 265;
CAPTURE(actual);
CAPTURE(expected);
CHECK(std::abs(actual-expected) < 0.01);
}
SECTION("Ice VI-liquid")
{
double actual = AS->melting_line(iT, iP, 1356.76e6);
double expected = 320;
CAPTURE(actual);
CAPTURE(expected);
CHECK(std::abs(actual-expected) < 1);
}
}
#endif

7
src/Tests/CoolProp-Tests.cpp
View File

@@ -778,8 +778,8 @@ TEST_CASE("Tests for solvers in P,Y flash using Water", "[flash],[PH],[PS],[PU]"
{
double Tc = Props1SI("Water","Tcrit");
double pc = Props1SI("Water","pcrit");
double p = pc*1.3;
double T = Tc*0.7;
double p = pc*2;
double T = Tc*0.5;
CAPTURE(T);
CAPTURE(p);
CHECK(ValidNumber(T));
@@ -793,9 +793,6 @@ TEST_CASE("Tests for solvers in P,Y flash using Water", "[flash],[PH],[PS],[PU]"
CHECK(ValidNumber(T2));
}
}
}
#endif

23
src/main.cxx
View File

@@ -18,7 +18,7 @@ using namespace CoolProp;
#endif
#include "SpeedTest.h"
#include "HumidAirProp.h"
//#include "CoolPropLib.h"
#include "CoolPropLib.h"
#include "crossplatform_shared_ptr.h"
@@ -400,23 +400,36 @@ int main()
}
#endif
#if 0
{
double Tc = Props1SI("Water","Tcrit");
double pc = Props1SI("Water","pcrit");
double p = pc*2;
double T = Tc*0.5;
char ykey[] = "H";
double y = PropsSI(ykey,"P",p,"T",T,"Water");
double TT = PropsSI("T","P",p,ykey,y,"Water");
int rr = 0;
}
#endif
#if 1
{
run_tests();
char c;
std::cin >> c;
}
#endif
#if 1
#if 0
{
char ykey[] = "H";
double Ts, y, T2, dT = -1;
double dd = PropsSI("T","P",101325,"T",114.357,"n-Propane");
shared_ptr<AbstractState> AS(AbstractState::factory("HEOS","water"));
double ptt = AS->melting_line(iP, iT, 273.159);
double ptt = AS->melting_line(iT, iP, 138.268e6);
Ts = PropsSI("T","P",101325,"Q",0,"n-Propane");
Ts = PropsSI("H","T",841.225,"P",2.86832e+007,"Water");
std::cout << get_global_param_string("errstring");
y = PropsSI(ykey,"T",Ts+dT,"P",101325,"n-Propane");
T2 = PropsSI("T",ykey,y,"P",101325,"n-Propane");