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

saturation p and saturation T work with Q = 0 for mixtures with phase envelopes

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-09-13 23:25:06 +02:00
parent cabe7a8a22
commit 260512b90f
2 changed files with 76 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -257,6 +257,17 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
PhaseEnvelopeData &env = HEOS.PhaseEnvelope;
std::size_t quality;
if (std::abs(HEOS._Q) < 100*DBL_EPSILON){
quality = 0;
}
else if (std::abs(HEOS._Q - 1) < 100*DBL_EPSILON){
quality = 1;
}
else{
throw ValueError("Quality is neither 0 or 1");
}
// 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){
@@ -265,25 +276,40 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
}
}
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);
if (quality == 1){
IO.bubble_point = false;
IO.y = HEOS.get_mole_fractions();
IO.x.resize(IO.y.size());
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);
}
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{
IO.bubble_point = true;
IO.x = HEOS.get_mole_fractions();
IO.y.resize(IO.x.size());
for (std::size_t i = 0; i < IO.y.size()-1; ++i) // First N-1 elements
{
IO.y[i] = CubicInterp(env.rhomolar_vap, env.x[i], imax-1, imax, imax+1, imax+2, IO.rhomolar_vap);
}
IO.y[IO.y.size()-1] = 1 - std::accumulate(IO.y.begin(), IO.y.end()-1, 0.0);
NR.call(HEOS, IO.x, IO.y, IO);
}
}
else if (solutions.size() == 0){
throw ValueError("No solution was found in PQ_flash");
@@ -404,6 +430,17 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
PhaseEnvelopeData &env = HEOS.PhaseEnvelope;
std::size_t quality;
if (std::abs(HEOS._Q) < 100*DBL_EPSILON){
quality = 0;
}
else if (std::abs(HEOS._Q - 1) < 100*DBL_EPSILON){
quality = 1;
}
else{
throw ValueError("Quality is neither 0 or 1");
}
// 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){
@@ -412,17 +449,21 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
}
}
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
if (quality == 1){
IO.bubble_point = false;
}
else{
IO.bubble_point = true;
}
IO.p = HEOS._p;
IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::P_IMPOSED;
// p -> rhomolar_vap
IO.rhomolar_vap = CubicInterp(env.p, env.rhomolar_vap, imax-1, imax, imax+1, imax+2, static_cast<long double>(IO.p));
IO.y = HEOS.get_mole_fractions(); // because Q = 1
IO.x = IO.y; // Just to give it good size
IO.x.resize(IO.y.size());
IO.T = CubicInterp(env.rhomolar_vap, env.T, 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

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

@@ -948,26 +948,31 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke
if (bubble_point){
for (unsigned int i = 0; i < N-1; ++i){
err_rel[i] = v[i]/y[i];
y[i] += v[i];
}
y[N-1] = 1 - std::accumulate(y.begin(), y.end()-1, 0.0);
}
else{
for (unsigned int i = 0; i < N-1; ++i){
err_rel[i] = v[i]/y[i];
err_rel[i] = v[i]/x[i];
x[i] += v[i];
}
x[N-1] = 1 - std::accumulate(x.begin(), x.end()-1, 0.0);
}
if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){
if (imposed_variable == newton_raphson_saturation_options::P_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){
else if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){
T += v[N-1]; err_rel[N-1] = v[N-1]/T;
rhomolar_liq += v[N]; err_rel[N] = v[N]/rhomolar_liq;
}
else{
throw ValueError("invalid imposed_variable");
}
//std::cout << format("\t%Lg ", this->error_rms) << T << " " << rhomolar_liq << " " << rhomolar_vap << " v " << vec_to_string(v, "%0.10Lg") << " x " << vec_to_string(x, "%0.10Lg") << " r " << vec_to_string(r, "%0.10Lg") << std::endl;
min_rel_change = min_abs_value(err_rel);
@@ -1003,12 +1008,15 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
if (bubble_point){
// Vapor is incipient phase, set its composition
rSatV.set_mole_fractions(y);
rSatL.set_mole_fractions(x);
}
else{
// Liquid is incipient phase, set its composition
rSatL.set_mole_fractions(x);
rSatV.set_mole_fractions(y);
}
std::vector<long double> &x = rSatL.get_mole_fractions();
std::vector<long double> &y = rSatV.get_mole_fractions();
if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){
rSatL.update(DmolarT_INPUTS, rhomolar_liq, T);
rSatV.update(DmolarT_INPUTS, rhomolar_vap, T);
@@ -1041,13 +1049,11 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
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
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
if (bubble_point){
J[i][j] = -MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatV, i, j, xN_flag);
}
}
else{
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
else{
J[i][j] = MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatL, i, j, xN_flag);
}
}
@@ -1066,7 +1072,7 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
J[N][N] = rSatL.first_partial_deriv(iP, iDmolar, iT);
}
else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED){
// Independent variables are N-1 mole fractions of incipient phase, T and rho'
// Independent variables are N-1 mole fractions of incipient phase and T
// For the residuals F_i (equality of fugacities)
for (std::size_t i = 0; i < N; ++i)
@@ -1076,13 +1082,11 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
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
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
if (bubble_point){
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
else{
J[i][j] = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
}
}
@@ -1090,7 +1094,7 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
}
}
else if (imposed_variable == newton_raphson_saturation_options::T_IMPOSED){
// Independent variables are N-1 mole fractions of incipient phase, p and rho'
// Independent variables are N-1 mole fractions of incipient phase and p
// For the residuals F_i (equality of fugacities)
for (std::size_t i = 0; i < N; ++i)
@@ -1100,13 +1104,11 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
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
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
if (bubble_point){
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
else{
J[i][j] = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
}
}