Implemented all phase functions for BICUBIC; something funny with TTSE still

closes #763
closes #789

Let me know if I missed something
This commit is contained in:
Ian Bell
2015-10-02 16:35:53 -06:00
parent 915416d47c
commit 80851f4c49
5 changed files with 171 additions and 34 deletions

View File

@@ -39,11 +39,20 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
std::size_t iL, iV, iclosest = 0;
CoolPropDbl hL = 0, hV = 0;
SimpleState closest_state;
if (
(is_mixture && PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, iHmolar, _hmolar, iclosest, closest_state))
||
(!is_mixture && pure_saturation.is_inside(iP, _p, iHmolar, _hmolar, iL, iV, hL, hV))
)
bool is_two_phase = false;
// Phase is imposed, use it
if (imposed_phase_index != iphase_not_imposed){
is_two_phase = (imposed_phase_index == iphase_twophase);
}
else{
if (is_mixture){
is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, iHmolar, _hmolar, iclosest, closest_state);
}
else{
is_two_phase = pure_saturation.is_inside(iP, _p, iHmolar, _hmolar, iL, iV, hL, hV);
}
}
if ( is_two_phase )
{
using_single_phase_table = false;
_Q = (static_cast<double>(_hmolar)-hL)/(hV-hL);
@@ -52,6 +61,7 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
}
else{
cached_saturation_iL = iL; cached_saturation_iV = iV;
_phase = iphase_twophase;
}
}
else{
@@ -68,6 +78,8 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
if (!cell.valid()){throw ValueError(format("Cell is invalid and has no good neighbors for hmolar = %g, p= %g",val1,val2));}
}
}
// Recalculate the phase
recalculate_singlephase_phase();
}
}
break;
@@ -92,11 +104,20 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
CoolPropDbl zL = 0, zV = 0;
std::size_t iclosest = 0;
SimpleState closest_state;
if (
(is_mixture && PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, otherkey, otherval, iclosest, closest_state))
||
(!is_mixture && pure_saturation.is_inside(iP, _p, otherkey, otherval, iL, iV, zL, zV))
){
bool is_two_phase = false;
// Phase is imposed, use it
if (imposed_phase_index != iphase_not_imposed){
is_two_phase = (imposed_phase_index == iphase_twophase);
}
else{
if (is_mixture){
is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, otherkey, otherval, iclosest, closest_state);
}
else{
is_two_phase = pure_saturation.is_inside(iP, _p, otherkey, otherval, iL, iV, zL, zV);
}
}
if ( is_two_phase ){
using_single_phase_table = false;
if (otherkey == iDmolar){
_Q = (1/otherval - 1/zL)/(1/zV - 1/zL);
@@ -110,6 +131,7 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
else{
cached_saturation_iL = iL; cached_saturation_iV = iV;
}
_phase = iphase_twophase;
}
else{
// Find and cache the indices i, j
@@ -127,6 +149,8 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
}
// Now find hmolar given P, X for X in Hmolar, Smolar, Umolar
invert_single_phase_x(single_phase_logph, dataset->coeffs_ph, otherkey, otherval, _p, cached_single_phase_i, cached_single_phase_j);
// Recalculate the phase
recalculate_singlephase_phase();
}
break;
}
@@ -142,7 +166,7 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
// Call again, but this time with molar units; S: [J/kg/K] * [kg/mol] -> [J/mol/K]
update(PSmolar_INPUTS, val1, val2*AS->molar_mass()); return;
}
case PT_INPUTS:{
case PT_INPUTS:{
_p = val1; _T = val2;
if (!single_phase_logpT.native_inputs_are_in_range(_T, _p)){
// Use the AbstractState instance
@@ -155,11 +179,20 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
std::size_t iL = 0, iV = 0, iclosest = 0;
CoolPropDbl TL = 0, TV = 0;
SimpleState closest_state;
if (
(is_mixture && PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, iT, _T, iclosest, closest_state))
||
(!is_mixture && pure_saturation.is_inside(iP, _p, iT, _T, iL, iV, TL, TV))
)
bool is_two_phase = false;
// Phase is imposed, use it
if (imposed_phase_index != iphase_not_imposed){
is_two_phase = (imposed_phase_index == iphase_twophase);
}
else{
if (is_mixture){
is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, iT, _T, iclosest, closest_state);
}
else{
is_two_phase = pure_saturation.is_inside(iP, _p, iT, _T, iL, iV, TL, TV);
}
}
if ( is_two_phase )
{
using_single_phase_table = false;
throw ValueError(format("P,T with TTSE cannot be two-phase for now"));
@@ -199,6 +232,8 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
}
}
}
// Recalculate the phase
recalculate_singlephase_phase();
}
}
break;
@@ -225,11 +260,20 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
CoolPropDbl zL = 0, zV = 0;
std::size_t iclosest = 0;
SimpleState closest_state;
if (
(is_mixture && PhaseEnvelopeRoutines::is_inside(phase_envelope, iT, _T, otherkey, otherval, iclosest, closest_state))
||
(!is_mixture && pure_saturation.is_inside(iT, _T, otherkey, otherval, iL, iV, zL, zV))
)
bool is_two_phase = false;
// Phase is imposed, use it
if (imposed_phase_index != iphase_not_imposed){
is_two_phase = (imposed_phase_index == iphase_twophase);
}
else{
if (is_mixture){
is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iT, _T, otherkey, otherval, iclosest, closest_state);
}
else{
is_two_phase = pure_saturation.is_inside(iT, _T, otherkey, otherval, iL, iV, zL, zV);
}
}
if ( is_two_phase )
{
using_single_phase_table = false;
if (otherkey == iDmolar){
@@ -262,6 +306,8 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
}
// Now find the y variable (Dmolar or Smolar in this case)
invert_single_phase_y(single_phase_logpT, dataset->coeffs_pT, otherkey, otherval, _T, cached_single_phase_i, cached_single_phase_j);
// Recalculate the phase
recalculate_singlephase_phase();
}
break;
}
@@ -286,7 +332,7 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
}
}
_T = _Q*TV + (1-_Q)*TL;
cached_saturation_iL = iL; cached_saturation_iV = iV;
cached_saturation_iL = iL; cached_saturation_iV = iV; _phase = iphase_twophase;
}
break;
}
@@ -312,7 +358,7 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
pure_saturation.is_inside(iT, _T, iQ, _Q, iL, iV, pL, pV);
}
_p = _Q*pV + (1-_Q)*pL;
cached_saturation_iL = iL; cached_saturation_iV = iV;
cached_saturation_iL = iL; cached_saturation_iV = iV; _phase = iphase_twophase;
}
break;
}

View File

@@ -63,11 +63,10 @@ A^{-1} = \left[ \begin{array}{*{16}c} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
typedef std::vector<std::vector<double> > mat;
class BicubicBackend : public TabularBackend
{
protected:
public:
/// Instantiator; base class loads or makes tables
BicubicBackend(shared_ptr<CoolProp::AbstractState> AS) : TabularBackend (AS){
BicubicBackend(shared_ptr<CoolProp::AbstractState> AS) : TabularBackend(AS){
imposed_phase_index = iphase_not_imposed;
// If a pure fluid, don't need to set fractions, go ahead and build
if (this->AS->get_mole_fractions().size() == 1){
check_tables();

View File

@@ -39,20 +39,37 @@ void CoolProp::TTSEBackend::update(CoolProp::input_pairs input_pair, double val1
using_single_phase_table = true; // Use the table !
std::size_t iL, iV;
CoolPropDbl hL = 0, hV = 0;
if (pure_saturation.is_inside(iP, _p, iHmolar, _hmolar, iL, iV, hL, hV)){
std::size_t iclosest = 0;
SimpleState closest_state;
bool is_two_phase = false;
// Phase is imposed, use it
if (imposed_phase_index != iphase_not_imposed){
is_two_phase = (imposed_phase_index == iphase_twophase);
}
else{
if (is_mixture){
is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, iHmolar, _hmolar, iclosest, closest_state);
}
else{
is_two_phase = pure_saturation.is_inside(iP, _p, iHmolar, _hmolar, iL, iV, hL, hV);
}
}
if ( is_two_phase ){
using_single_phase_table = false;
_Q = (static_cast<double>(_hmolar)-hL)/(hV-hL);
if(!is_in_closed_range(0.0,1.0,static_cast<double>(_Q))){
throw ValueError("vapor quality is not in (0,1)");
}
else{
cached_saturation_iL = iL; cached_saturation_iV = iV;
cached_saturation_iL = iL; cached_saturation_iV = iV; _phase = iphase_twophase;
}
}
else{
// Find and cache the indices i, j
selected_table = SELECTED_PH_TABLE;
single_phase_logph.find_native_nearest_good_neighbor(_hmolar, _p, cached_single_phase_i, cached_single_phase_j);
// Recalculate the phase
recalculate_singlephase_phase();
}
}
break;
@@ -75,7 +92,23 @@ void CoolProp::TTSEBackend::update(CoolProp::input_pairs input_pair, double val1
using_single_phase_table = true; // Use the table (or first guess is that it is single-phase)!
std::size_t iL, iV;
CoolPropDbl zL = 0, zV = 0;
if (pure_saturation.is_inside(iP, _p, otherkey, otherval, iL, iV, zL, zV)){
std::size_t iclosest = 0;
SimpleState closest_state;
bool is_two_phase = false;
// Phase is imposed, use it
if (imposed_phase_index != iphase_not_imposed)
{
is_two_phase = (imposed_phase_index == iphase_twophase);
}
else{
if (is_mixture){
is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iP, _p, otherkey, otherval, iclosest, closest_state);
}
else{
is_two_phase = pure_saturation.is_inside(iP, _p, otherkey, otherval, iL, iV, zL, zV);
}
}
if (is_two_phase){
using_single_phase_table = false;
if (otherkey == iDmolar){
_Q = (1/otherval - 1/zL)/(1/zV - 1/zL);
@@ -87,7 +120,7 @@ void CoolProp::TTSEBackend::update(CoolProp::input_pairs input_pair, double val1
throw ValueError("vapor quality is not in (0,1)");
}
else{
cached_saturation_iL = iL; cached_saturation_iV = iV;
cached_saturation_iL = iL; cached_saturation_iV = iV; _phase = iphase_twophase;
}
}
else{
@@ -96,6 +129,8 @@ void CoolProp::TTSEBackend::update(CoolProp::input_pairs input_pair, double val1
single_phase_logph.find_nearest_neighbor(iP, _p, otherkey, otherval, cached_single_phase_i, cached_single_phase_j);
// Now find hmolar
invert_single_phase_x(single_phase_logph, otherkey, otherval, _p, cached_single_phase_i, cached_single_phase_j);
// Recalculate the phase
recalculate_singlephase_phase();
}
break;
}
@@ -131,6 +166,8 @@ void CoolProp::TTSEBackend::update(CoolProp::input_pairs input_pair, double val1
// Find and cache the indices i, j
selected_table = SELECTED_PT_TABLE;
single_phase_logpT.find_native_nearest_neighbor(_T, _p, cached_single_phase_i, cached_single_phase_j);
// Recalculate the phase
recalculate_singlephase_phase();
}
}
break;
@@ -155,7 +192,22 @@ void CoolProp::TTSEBackend::update(CoolProp::input_pairs input_pair, double val1
using_single_phase_table = true; // Use the table (or first guess is that it is single-phase)!
std::size_t iL, iV;
CoolPropDbl zL = 0, zV = 0;
if (pure_saturation.is_inside(iT, _T, otherkey, otherval, iL, iV, zL, zV)){
std::size_t iclosest = 0;
SimpleState closest_state;
bool is_two_phase = false;
// Phase is imposed, use it
if (imposed_phase_index != iphase_not_imposed){
is_two_phase = (imposed_phase_index == iphase_twophase);
}
else{
if (is_mixture){
is_two_phase = PhaseEnvelopeRoutines::is_inside(phase_envelope, iT, _T, otherkey, otherval, iclosest, closest_state);
}
else{
is_two_phase = pure_saturation.is_inside(iT, _T, otherkey, otherval, iL, iV, zL, zV);
}
}
if ( is_two_phase ){
using_single_phase_table = false;
if (otherkey == iDmolar){
_Q = (1/otherval - 1/zL)/(1/zV - 1/zL);
@@ -167,7 +219,7 @@ void CoolProp::TTSEBackend::update(CoolProp::input_pairs input_pair, double val1
throw ValueError("vapor quality is not in (0,1)");
}
else{
cached_saturation_iL = iL; cached_saturation_iV = iV;
cached_saturation_iL = iL; cached_saturation_iV = iV; _phase = iphase_twophase;
}
}
else{
@@ -176,6 +228,8 @@ void CoolProp::TTSEBackend::update(CoolProp::input_pairs input_pair, double val1
single_phase_logpT.find_nearest_neighbor(iT, _T, otherkey, otherval, cached_single_phase_i, cached_single_phase_j);
// Now find hmolar
invert_single_phase_y(single_phase_logpT, otherkey, otherval, _T, cached_single_phase_i, cached_single_phase_j);
// Recalculate the phase
recalculate_singlephase_phase();
}
break;
}
@@ -200,7 +254,7 @@ void CoolProp::TTSEBackend::update(CoolProp::input_pairs input_pair, double val1
}
}
_T = _Q*TV + (1-_Q)*TL;
cached_saturation_iL = iL; cached_saturation_iV = iV;
cached_saturation_iL = iL; cached_saturation_iV = iV; _phase = iphase_twophase;
}
break;
}
@@ -226,7 +280,7 @@ void CoolProp::TTSEBackend::update(CoolProp::input_pairs input_pair, double val1
pure_saturation.is_inside(iT, _T, iQ, _Q, iL, iV, pL, pV);
}
_p = _Q*pV + (1-_Q)*pL;
cached_saturation_iL = iL; cached_saturation_iV = iV;
cached_saturation_iL = iL; cached_saturation_iV = iV; _phase = iphase_twophase;
}
break;
}

View File

@@ -12,6 +12,7 @@ class TTSEBackend : public TabularBackend
std::string backend_name(void){return "TTSEBackend";}
/// Instantiator; base class loads or makes tables
TTSEBackend(shared_ptr<CoolProp::AbstractState> AS) : TabularBackend (AS) {
imposed_phase_index = iphase_not_imposed;
// If a pure fluid, don't need to set fractions, go ahead and build
if (this->AS->get_mole_fractions().size() == 1){
check_tables();

View File

@@ -706,6 +706,7 @@ public:
class TabularBackend : public AbstractState
{
protected:
phases imposed_phase_index;
bool tables_loaded, using_single_phase_table, is_mixture;
enum selected_table_options{SELECTED_NO_TABLE=0, SELECTED_PH_TABLE, SELECTED_PT_TABLE};
selected_table_options selected_table;
@@ -762,6 +763,42 @@ class TabularBackend : public AbstractState
}
TabularDataSet * dataset;
void recalculate_singlephase_phase()
{
if (p() > p_critical()){
if (T() > T_critical()){
_phase = iphase_supercritical;
}
else{
_phase = iphase_supercritical_liquid;
}
}
else{
if (T() > T_critical()){
_phase = iphase_supercritical_gas;
}
else{
// Liquid or vapor
if (rhomolar() > rhomolar_critical()){
_phase = iphase_liquid;
}
else{
_phase = iphase_gas;
}
}
}
}
/** \brief Specify the phase - this phase will always be used in calculations
*
* @param phase_index The index from CoolProp::phases
*/
void calc_specify_phase(phases phase_index){ imposed_phase_index = phase_index; };
/**\brief Unspecify the phase - the phase is no longer imposed, different solvers can do as they like
*/
void calc_unspecify_phase(){ imposed_phase_index = iphase_not_imposed; };
phases calc_phase(void){ return _phase; }
CoolPropDbl calc_T_critical(void){return this->AS->T_critical();};
CoolPropDbl calc_Ttriple(void){return this->AS->Ttriple();};
CoolPropDbl calc_p_triple(void){return this->AS->p_triple();};