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

QT calculations using phase envelope work for Q = 1

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-09-13 21:42:31 +02:00
parent 8c843bbf18
commit cabe7a8a22
3 changed files with 94 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -251,22 +251,70 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
}
else
{
// Set some input options
SaturationSolvers::mixture_VLE_IO options;
options.sstype = SaturationSolvers::imposed_T;
options.Nstep_max = 20;
if(HEOS.PhaseEnvelope.built){
// Find the intersections in the phase envelope
std::vector< std::pair<std::size_t, std::size_t> > intersections = PhaseEnvelopeRoutines::find_intersections(HEOS, iT, HEOS._T);
PhaseEnvelopeData &env = HEOS.PhaseEnvelope;
// Find the correct solution
std::vector<std::size_t> solutions;
for (std::vector< std::pair<std::size_t, std::size_t> >::iterator it = intersections.begin(); it != intersections.end(); ++it){
if (std::abs(env.Q[it->first] - HEOS._Q) < 10*DBL_EPSILON && std::abs(env.Q[it->second] - HEOS._Q) < 10*DBL_EPSILON ){
solutions.push_back(it->first);
}
}
if (solutions.size() == 1){
if (std::abs(HEOS._Q-1) > 1e-6){ throw NotImplementedError("Q = 1 only for now"); }
std::size_t &imax = solutions[0];
SaturationSolvers::newton_raphson_saturation NR;
SaturationSolvers::newton_raphson_saturation_options IO;
IO.bubble_point = false; // because Q = 1
IO.T = HEOS._T;
IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::T_IMPOSED;
// p -> rhomolar_vap
IO.rhomolar_vap = CubicInterp(env.T, env.rhomolar_vap, imax-1, imax, imax+1, imax+2, static_cast<long double>(IO.T));
IO.y = HEOS.get_mole_fractions(); // because Q = 1
IO.x = IO.y; // Just to give it good size
IO.p = CubicInterp(env.rhomolar_vap, env.p, imax-1, imax, imax+1, imax+2, IO.rhomolar_vap);
IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, imax-1, imax, imax+1, imax+2, IO.rhomolar_vap);
for (std::size_t i = 0; i < IO.x.size()-1; ++i) // First N-1 elements
{
IO.x[i] = CubicInterp(env.rhomolar_vap, env.x[i], imax-1, imax, imax+1, imax+2, IO.rhomolar_vap);
}
IO.x[IO.x.size()-1] = 1 - std::accumulate(IO.x.begin(), IO.x.end()-1, 0.0);
NR.call(HEOS, IO.y, IO.x, IO);
}
else if (solutions.size() == 0){
throw ValueError("No solution was found in PQ_flash");
}
else{
throw ValueError("More than 1 solution was ");
}
// Load the outputs
HEOS._phase = iphase_twophase;
HEOS._p = HEOS.SatV->p();
HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar());
HEOS._T = HEOS.SatL->T();
}
else{
// Set some input options
SaturationSolvers::mixture_VLE_IO options;
options.sstype = SaturationSolvers::imposed_T;
options.Nstep_max = 20;
// Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted
long double pguess = SaturationSolvers::saturation_preconditioner(HEOS, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions);
// Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted
long double pguess = SaturationSolvers::saturation_preconditioner(HEOS, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions);
// Use Wilson iteration to obtain updated guess for pressure
pguess = SaturationSolvers::saturation_Wilson(HEOS, HEOS._Q, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions, pguess);
// Use Wilson iteration to obtain updated guess for pressure
pguess = SaturationSolvers::saturation_Wilson(HEOS, HEOS._Q, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions, pguess);
// Actually call the successive substitution solver
SaturationSolvers::successive_substitution(HEOS, HEOS._Q, HEOS._T, pguess, HEOS.mole_fractions, HEOS.K, options);
// Actually call the successive substitution solver
SaturationSolvers::successive_substitution(HEOS, HEOS._Q, HEOS._T, pguess, HEOS.mole_fractions, HEOS.K, options);
HEOS._p = options.p;
HEOS._rhomolar = 1/(HEOS._Q/options.rhomolar_vap+(1-HEOS._Q)/options.rhomolar_liq);
HEOS._p = options.p;
HEOS._rhomolar = 1/(HEOS._Q/options.rhomolar_vap+(1-HEOS._Q)/options.rhomolar_liq);
}
}
}
void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)

37
src/Backends/Helmholtz/VLERoutines.cpp
View File

@@ -805,7 +805,7 @@ void SaturationSolvers::newton_raphson_saturation::resize(unsigned int N)
}
err_rel.resize(N+1);
}
else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED){
else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::T_IMPOSED){
r.resize(N);
negative_r.resize(N);
J.resize(N, std::vector<long double>(N, 0));
@@ -959,7 +959,12 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke
}
x[N-1] = 1 - std::accumulate(x.begin(), x.end()-1, 0.0);
}
T += v[N-1]; err_rel[N-1] = v[N-1]/T;
if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){
T += v[N-1]; err_rel[N-1] = v[N-1]/T;
}
else if (imposed_variable == newton_raphson_saturation_options::T_IMPOSED){
p += v[N-1]; err_rel[N-1] = v[N-1]/p;
}
if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){
rhomolar_liq += v[N]; err_rel[N] = v[N]/rhomolar_liq;
}
@@ -1008,9 +1013,9 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
rSatL.update(DmolarT_INPUTS, rhomolar_liq, T);
rSatV.update(DmolarT_INPUTS, rhomolar_vap, T);
}
else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED){
else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::T_IMPOSED){
rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar();
rSatV.update_TP_guessrho(T, p, rhomolar_vap);
rSatV.update_TP_guessrho(T, p, rhomolar_vap); rhomolar_vap = rSatV.rhomolar();
}
else{
throw ValueError("imposed variable not set for NR VLE");
@@ -1084,6 +1089,30 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
J[i][N-1] = MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatV, i, xN_flag);
}
}
else if (imposed_variable == newton_raphson_saturation_options::T_IMPOSED){
// Independent variables are N-1 mole fractions of incipient phase, p and rho'
// For the residuals F_i (equality of fugacities)
for (std::size_t i = 0; i < N; ++i)
{
// Equate the liquid and vapor fugacities
long double ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
long double ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
r[i] = ln_f_liq - ln_f_vap;
if (bubble_point){
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
J[i][j] = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag);
}
}
else{
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
J[i][j] = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
}
}
J[i][N-1] = MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatV, i, xN_flag);
}
}
else{
throw ValueError();
}

2
src/Backends/Helmholtz/VLERoutines.h
View File

@@ -201,7 +201,7 @@ namespace SaturationSolvers
};
struct newton_raphson_saturation_options{
enum imposed_variable_options {NO_VARIABLE_IMPOSED = 0, P_IMPOSED, RHOV_IMPOSED};
enum imposed_variable_options {NO_VARIABLE_IMPOSED = 0, P_IMPOSED, RHOV_IMPOSED, T_IMPOSED};
int Nstep_max;
bool bubble_point;
std::size_t Nsteps;