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

Add (passing) tests for stability code, tests for critical point calcs for paper

Browse Source
This commit is contained in:
Ian Bell
2016-08-23 20:57:43 -06:00
parent 8f0aa03650
commit 05947657b1
2 changed files with 77 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -7,12 +7,11 @@
#if defined(ENABLE_CATCH)
#include "catch.hpp"
#include "Backends/Cubics/CubicBackend.h"
#endif
namespace CoolProp{
void FlashRoutines::PT_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.PhaseEnvelope.built){
@@ -1888,6 +1887,8 @@ void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS)
#if defined(ENABLE_CATCH)
TEST_CASE("PD with T very large should yield error","[PDflash]")
{
shared_ptr<HelmholtzEOSBackend> HEOS(new HelmholtzEOSBackend("R134a"));
@@ -1898,7 +1899,7 @@ TEST_CASE("PD with T very large should yield error","[PDflash]")
TEST_CASE("Stability testing","[stability]")
{
shared_ptr<HelmholtzEOSMixtureBackend> HEOS(new HelmholtzEOSMixtureBackend(strsplit("Methane&Ethane&n-Propane&n-Butane",'&')));
shared_ptr<HelmholtzEOSMixtureBackend> HEOS(new HelmholtzEOSMixtureBackend(strsplit("n-Propane&n-Butane&n-Pentane&n-Hexane",'&')));
std::vector<double> z(4); z[0] = 0.1; z[1] = 0.2; z[2] = 0.3; z[3] = 0.4;
HEOS->set_mole_fractions(z);
@@ -1910,13 +1911,21 @@ TEST_CASE("Stability testing","[stability]")
SECTION("Liquid (feed is stable)"){
StabilityRoutines::StabilityEvaluationClass stability_tester(*HEOS);
stability_tester.set_TP(TL-1, 101325);
CHECK(stability_tester.is_stable() == true);
for (double T = TL-1; T >= 100; T -= 1)
{
stability_tester.set_TP(T, 101325);
CAPTURE(T);
CHECK_NOTHROW(stability_tester.is_stable());
}
}
SECTION("Vapor (feed is stable)"){
StabilityRoutines::StabilityEvaluationClass stability_tester(*HEOS);
stability_tester.set_TP(TV+1, 101325);
CHECK(stability_tester.is_stable() == true);
for (double T = TV+1; T <= 500; T += 1)
{
stability_tester.set_TP(T, 101325);
CAPTURE(T);
CHECK_NOTHROW(stability_tester.is_stable());
}
}
SECTION("Two-phase (feed is unstable)"){
StabilityRoutines::StabilityEvaluationClass stability_tester(*HEOS);
@@ -1925,7 +1934,52 @@ TEST_CASE("Stability testing","[stability]")
}
}
TEST_CASE("Test critical points for methane + H2S","[critical_points]")
{
shared_ptr<HelmholtzEOSMixtureBackend> HEOS(new HelmholtzEOSMixtureBackend(strsplit("Methane&H2S",'&')));
double zz[] = {0.998, 0.97, 0.9475, 0.94, 0.93, 0.86, 0.85, 0.84, 0.75, 0.53, 0.52, 0.51, 0.49, 0.36, 0.24, 0.229, 0.09};
int Npts[] = {2,2,2,2,2,2,2,2,0,0,2,2,2,1,1,1,1}; // Number of critical points that should be found
int imax = sizeof(zz)/sizeof(double);
for (int i = 0; i < imax; ++i){
double z0 = zz[i];
std::vector<double> z(2); z[0] = z0; z[1] = 1-z0;
HEOS->set_mole_fractions(z);
CAPTURE(z0);
std::vector<CriticalState> pts = HEOS->all_critical_points();
CHECK(pts.size() == Npts[i]);
}
}
TEST_CASE("Test critical points for nitrogen + ethane with HEOS","[critical_points]")
{
shared_ptr<HelmholtzEOSMixtureBackend> HEOS(new HelmholtzEOSMixtureBackend(strsplit("Nitrogen&Ethane",'&')));
std::vector<double> zz = linspace(0.001, 0.999, 21);
for (int i = 0; i < static_cast<std::size_t>(zz.size()); ++i){
double z0 = zz[i];
std::vector<double> z(2); z[0] = z0; z[1] = 1-z0;
HEOS->set_mole_fractions(z);
CAPTURE(z0);
std::vector<CriticalState> pts;
CHECK_NOTHROW(pts = HEOS->all_critical_points(););
}
}
TEST_CASE("Test critical points for nitrogen + ethane with PR","[critical_points]")
{
shared_ptr<PengRobinsonBackend> HEOS(new PengRobinsonBackend(strsplit("Nitrogen&Ethane",'&')));
HEOS->set_binary_interaction_double(0, 1, "kij", 0.0407); // Ramırez-Jimenez et al.
std::vector<double> zz = linspace(0.001, 0.999, 21);
for (int i = 0; i < static_cast<std::size_t>(zz.size()); ++i){
double z0 = zz[i];
std::vector<double> z(2); z[0] = z0; z[1] = 1-z0;
HEOS->set_mole_fractions(z);
CAPTURE(z0);
std::vector<CriticalState> pts;
CHECK_NOTHROW(pts = HEOS->all_critical_points(););
}
}
#endif

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

@@ -2117,14 +2117,21 @@ HelmholtzEOSBackend::StationaryPointReturnFlag HelmholtzEOSMixtureBackend::solve
}
}catch(std::exception &e){ if (get_debug_level() > 5) {std::cout << e.what() << std::endl; }; light = -1;}
try{
heavy = Halley(resid, rhomax, 1e-8, 100);
double d2pdrho2__constT = resid.deriv(heavy);
if (d2pdrho2__constT < 0){
// Not possible since curvature should be positive
throw CoolProp::ValueError("curvature cannot be negative");
}
}catch(std::exception &e){ if (get_debug_level() > 5) {std::cout << e.what() << std::endl; }; heavy = -1;}
// First try a "normal" calculation of the stationary point on the liquid side
for (double omega = 1.0; omega > 0; omega -= 0.4){
try{
resid.options.add_number("omega", omega);
heavy = Halley(resid, rhomax, 1e-8, 100);
double d2pdrho2__constT = resid.deriv(heavy);
if (d2pdrho2__constT < 0){
// Not possible since curvature should be positive
throw CoolProp::ValueError("curvature cannot be negative");
}
break; // Jump out, we got a good solution
}catch(std::exception &e){ if (get_debug_level() > 5) {std::cout << e.what() << std::endl; }; heavy = -1;}
}
if (light > 0 && heavy > 0){
// Found two stationary points, done!
return TWO_STATIONARY_POINTS_FOUND;
@@ -2143,7 +2150,7 @@ HelmholtzEOSBackend::StationaryPointReturnFlag HelmholtzEOSMixtureBackend::solve
}
}
else{
throw CoolProp::ValueError("one stationary points -- ugh?");
return ONE_STATIONARY_POINT_FOUND;
}
}
// Define the residual to be driven to zero