github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-01-22 20:38:01 -05:00
Code Issues Packages Projects Releases Wiki Activity

Added solver for mixtures with vapor quality set and one of P or T. Q no longer needs to be 0 or 1

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-09-14 22:11:59 +02:00
parent 44f1c52a61
commit 20e55df16e
4 changed files with 424 additions and 140 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -252,76 +252,7 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
else
{
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;
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){
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){
std::size_t &imax = solutions[0];
SaturationSolvers::newton_raphson_saturation NR;
SaturationSolvers::newton_raphson_saturation_options IO;
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.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);
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);
}
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");
}
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();
PT_Q_flash_mixtures(HEOS, iT, HEOS._T);
}
else{
// Set some input options
@@ -341,6 +272,11 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
HEOS._p = options.p;
HEOS._rhomolar = 1/(HEOS._Q/options.rhomolar_vap+(1-HEOS._Q)/options.rhomolar_liq);
}
// 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();
}
}
void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
@@ -425,65 +361,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
else
{
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, iP, HEOS._p);
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){
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){
std::size_t &imax = solutions[0];
SaturationSolvers::newton_raphson_saturation NR;
SaturationSolvers::newton_raphson_saturation_options IO;
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.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
{
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();
PT_Q_flash_mixtures(HEOS, iP, HEOS._p);
}
else{
// Set some imput options
@@ -499,13 +377,189 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
// Actually call the successive substitution solver
SaturationSolvers::successive_substitution(HEOS, HEOS._Q, Tguess, HEOS._p, HEOS.mole_fractions, HEOS.K, io);
// 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();
}
// 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();
}
}
void FlashRoutines::PT_Q_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS, parameters other, long double value)
{
// Find the intersections in the phase envelope
std::vector< std::pair<std::size_t, std::size_t> > intersections = PhaseEnvelopeRoutines::find_intersections(HEOS, other, value);
PhaseEnvelopeData &env = HEOS.PhaseEnvelope;
enum quality_options{SATURATED_LIQUID, SATURATED_VAPOR, TWO_PHASE};
quality_options quality;
if (std::abs(HEOS._Q) < 100*DBL_EPSILON){
quality = SATURATED_LIQUID;
}
else if (std::abs(HEOS._Q - 1) < 100*DBL_EPSILON){
quality = SATURATED_VAPOR;
}
else if (HEOS._Q > 0 && HEOS._Q < 1){
quality = TWO_PHASE;
}
else{
throw ValueError("Quality is not within 0 and 1");
}
if (quality == SATURATED_LIQUID || quality == SATURATED_VAPOR)
{
// *********************************************************
// Bubble- or dew-point calculation
// *********************************************************
// 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){
std::size_t &imax = solutions[0];
SaturationSolvers::newton_raphson_saturation NR;
SaturationSolvers::newton_raphson_saturation_options IO;
if (other == iP){
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));
}
else if (other == iT){
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));
}
if (quality == SATURATED_VAPOR){
IO.bubble_point = false;
IO.y = HEOS.get_mole_fractions(); // Because Q = 1
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);
}
else{
IO.bubble_point = true;
IO.x = HEOS.get_mole_fractions(); // Because Q = 0
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");
}
else{
throw ValueError("More than 1 solution was found in PQ_flash");
}
}
else{
// *********************************************************
// Two-phase calculation for given vapor quality
// *********************************************************
// Find the correct solution
std::vector<std::size_t> liquid_solutions, vapor_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] - 0) < 10*DBL_EPSILON && std::abs(env.Q[it->second] - 0) < 10*DBL_EPSILON ){
liquid_solutions.push_back(it->first);
}
if (std::abs(env.Q[it->first] - 1) < 10*DBL_EPSILON && std::abs(env.Q[it->second] - 1) < 10*DBL_EPSILON ){
vapor_solutions.push_back(it->first);
}
}
if (liquid_solutions.size() != 1 || vapor_solutions.size() != 1){
throw ValueError(format("Number liquid solutions [%d] or vapor solutions [%d] != 1", liquid_solutions.size(), vapor_solutions.size() ));
}
std::size_t iliq = liquid_solutions[0], ivap = vapor_solutions[0];
SaturationSolvers::newton_raphson_twophase NR;
SaturationSolvers::newton_raphson_twophase_options IO;
IO.beta = HEOS._Q;
long double rhomolar_vap_sat_vap, T_sat_vap, rhomolar_liq_sat_vap, rhomolar_liq_sat_liq, T_sat_liq, rhomolar_vap_sat_liq, p_sat_liq, p_sat_vap;
if (other == iP){
IO.p = HEOS._p;
p_sat_liq = IO.p; p_sat_vap = IO.p;
IO.imposed_variable = SaturationSolvers::newton_raphson_twophase_options::P_IMPOSED;
// Calculate the interpolated values for beta = 0 and beta = 1
rhomolar_vap_sat_vap = CubicInterp(env.p, env.rhomolar_vap, ivap-1, ivap, ivap+1, ivap+2, static_cast<long double>(IO.p));
T_sat_vap = CubicInterp(env.rhomolar_vap, env.T, ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap);
rhomolar_liq_sat_vap = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap);
// Phase inversion for liquid solution (liquid is vapor and vice versa)
rhomolar_liq_sat_liq = CubicInterp(env.p, env.rhomolar_vap, iliq-1, iliq, iliq+1, iliq+2, static_cast<long double>(IO.p));
T_sat_liq = CubicInterp(env.rhomolar_vap, env.T, iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq);
rhomolar_vap_sat_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq);
}
else if (other == iT){
IO.T = HEOS._T;
T_sat_liq = IO.T; T_sat_vap = IO.T;
IO.imposed_variable = SaturationSolvers::newton_raphson_twophase_options::T_IMPOSED;
// Calculate the interpolated values for beta = 0 and beta = 1
rhomolar_vap_sat_vap = CubicInterp(env.T, env.rhomolar_vap, ivap-1, ivap, ivap+1, ivap+2, static_cast<long double>(IO.T));
p_sat_vap = CubicInterp(env.rhomolar_vap, env.p, ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap);
rhomolar_liq_sat_vap = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap);
// Phase inversion for liquid solution (liquid is vapor and vice versa)
rhomolar_liq_sat_liq = CubicInterp(env.T, env.rhomolar_vap, iliq-1, iliq, iliq+1, iliq+2, static_cast<long double>(IO.T));
p_sat_liq = CubicInterp(env.rhomolar_vap, env.p, iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq);
rhomolar_vap_sat_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq);
}
else{
throw ValueError();
}
// Weight the guesses by the vapor mole fraction
IO.rhomolar_vap = IO.beta*rhomolar_vap_sat_vap + (1-IO.beta)*rhomolar_vap_sat_liq;
IO.rhomolar_liq = IO.beta*rhomolar_liq_sat_vap + (1-IO.beta)*rhomolar_liq_sat_liq;
IO.T = IO.beta*T_sat_vap + (1-IO.beta)*T_sat_liq;
IO.p = IO.beta*p_sat_vap + (1-IO.beta)*p_sat_liq;
IO.z = HEOS.get_mole_fractions();
IO.x.resize(IO.z.size());
IO.y.resize(IO.z.size());
for (std::size_t i = 0; i < IO.x.size()-1; ++i) // First N-1 elements
{
long double x_sat_vap = CubicInterp(env.rhomolar_vap, env.x[i], ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap);
long double y_sat_vap = CubicInterp(env.rhomolar_vap, env.y[i], ivap-1, ivap, ivap+1, ivap+2, rhomolar_vap_sat_vap);
long double x_sat_liq = CubicInterp(env.rhomolar_vap, env.y[i], iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq);
long double y_sat_liq = CubicInterp(env.rhomolar_vap, env.x[i], iliq-1, iliq, iliq+1, iliq+2, rhomolar_liq_sat_liq);
IO.x[i] = IO.beta*x_sat_vap + (1-IO.beta)*x_sat_liq;
IO.y[i] = IO.beta*y_sat_vap + (1-IO.beta)*y_sat_liq;
}
IO.x[IO.x.size()-1] = 1 - std::accumulate(IO.x.begin(), IO.x.end()-1, 0.0);
IO.y[IO.y.size()-1] = 1 - std::accumulate(IO.y.begin(), IO.y.end()-1, 0.0);
std::vector<long double> &XX = IO.x;
std::vector<long double> &YY = IO.y;
NR.call(HEOS, IO);
}
}
// D given and one of P,H,S,U

6
src/Backends/Helmholtz/FlashRoutines.h
View File

@@ -36,6 +36,12 @@ public:
/// @param HEOS The HelmholtzEOSMixtureBackend to be used
static void QT_flash(HelmholtzEOSMixtureBackend &HEOS);
/// Flash for mixture given temperature or pressure and (molar) quality
/// @param HEOS The HelmholtzEOSMixtureBackend to be used
/// @param other The parameter that is imposed, either iT or iP
/// @param value The value for the imposed parameter
static void PT_Q_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS, parameters other, long double value);
/// Flash for given pressure and temperature
/// @param HEOS The HelmholtzEOSMixtureBackend to be used
static void PT_flash(HelmholtzEOSMixtureBackend &HEOS);

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

@@ -1141,4 +1141,169 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
dPsat_dTsat = -dQ_dTsat/dQ_dPsat;
}
void SaturationSolvers::newton_raphson_twophase::call(HelmholtzEOSMixtureBackend &HEOS, newton_raphson_twophase_options &IO)
{
int iter = 0;
if (get_debug_level() > 9){std::cout << " NRsat::call: p" << IO.p << " T" << IO.T << " dl" << IO.rhomolar_liq << " dv" << IO.rhomolar_vap << std::endl;}
// Reset all the variables and resize
pre_call();
rhomolar_liq = IO.rhomolar_liq;
rhomolar_vap = IO.rhomolar_vap;
T = IO.T;
p = IO.p;
imposed_variable = IO.imposed_variable;
x = IO.x;
y = IO.y;
z = IO.z;
beta = IO.beta;
this->N = z.size();
x.resize(N);
y.resize(N);
r.resize(2*N-1);
negative_r.resize(2*N-1);
J.resize(2*N-1, std::vector<long double>(2*N-1, 0));
err_rel.resize(2*N-1);
// Hold a pointer to the backend
this->HEOS = &HEOS;
do
{
// Build the Jacobian and residual vectors
build_arrays();
// Solve for the step; v is the step with the contents
// [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)]
// Uncomment to see Jacobian and residual at every step
// std::cout << vec_to_string(J, "%0.12Lg") << std::endl;
// std::cout << vec_to_string(negative_r, "%0.12Lg") << std::endl;
std::vector<long double> v = linsolve(J, negative_r);
for (unsigned int i = 0; i < N-1; ++i){
err_rel[i] = v[i]/x[i];
x[i] += v[i];
err_rel[i+(N-1)] = v[i+(N-1)]/y[i];
y[i] += v[i+(N-1)];
}
x[N-1] = 1 - std::accumulate(x.begin(), x.end()-1, 0.0);
y[N-1] = 1 - std::accumulate(y.begin(), y.end()-1, 0.0);
if (imposed_variable == newton_raphson_twophase_options::P_IMPOSED){
T += v[2*N-2]; err_rel[2*N-2] = v[2*N-2]/T;
}
else if (imposed_variable == newton_raphson_twophase_options::T_IMPOSED){
p += v[2*N-2]; err_rel[2*N-2] = v[2*N-2]/p;
}
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);
iter++;
if (iter == IO.Nstep_max){
throw ValueError(format("newton_raphson_saturation::call reached max number of iterations [%d]",IO.Nstep_max));
}
}
while(this->error_rms > 1e-9 && min_rel_change > 1000*DBL_EPSILON && iter < IO.Nstep_max);
IO.Nsteps = iter;
IO.p = p;
IO.x = x; // Mole fractions in liquid
IO.y = y; // Mole fractions in vapor
IO.T = T;
IO.rhomolar_liq = rhomolar_liq;
IO.rhomolar_vap = rhomolar_vap;
IO.hmolar_liq = HEOS.SatL.get()->hmolar();
IO.hmolar_vap = HEOS.SatV.get()->hmolar();
IO.smolar_liq = HEOS.SatL.get()->smolar();
IO.smolar_vap = HEOS.SatV.get()->smolar();
}
void SaturationSolvers::newton_raphson_twophase::build_arrays()
{
// References to the classes for concision
HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
// Step 0:
// -------
// Set mole fractions
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();
rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar();
rSatV.update_TP_guessrho(T, p, rhomolar_vap); rhomolar_vap = rSatV.rhomolar();
// For diagnostic purposes calculate the pressures (no derivatives are evaluated)
long double p_liq = rSatL.p();
long double p_vap = rSatV.p();
p = 0.5*(p_liq + p_vap);
// Step 2:
// -------
// Build the residual vector and the Jacobian matrix
x_N_dependency_flag xN_flag = XN_DEPENDENT;
// Form of residuals do not depend on which variable is imposed
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; // N of these
if (i != N-1){
// Equate the specified vapor mole fraction and that given defined by the ith component
r[i+N] = (z[i]-x[i])/(y[i]-x[i]) - beta; // N-1 of these
}
}
// First part of derivatives with respect to ln f_i
for (std::size_t i = 0; i < N; ++i)
{
for (std::size_t j = 0; j < N-1; ++j){
J[i][j] = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
J[i][j+N-1] = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag);
}
// Last derivative with respect to either T or p depending on what is imposed
if (imposed_variable == newton_raphson_twophase_options::P_IMPOSED){
J[i][2*N-2] = 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_twophase_options::T_IMPOSED){
J[i][2*N-2] = 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();
}
}
// Derivatives with respect to the vapor mole fractions residual
for (std::size_t i = 0; i < N-1; ++i)
{
std::size_t k = i + N; // N ln f_i residuals
J[k][i] = (z[i]-y[i])/pow(y[i]-x[i], 2);
J[k][i+(N-1)] = -(z[i]-x[i])/pow(y[i]-x[i], 2);
}
// Flip all the signs of the entries in the residual vector since we are solving Jv = -r, not Jv=r
// Also calculate the rms error of the residual vector at this step
error_rms = 0;
for (unsigned int i = 0; i < J.size(); ++i)
{
negative_r[i] = -r[i];
error_rms += r[i]*r[i]; // Sum the squares
}
error_rms = sqrt(error_rms); // Square-root (The R in RMS)
}
} /* namespace CoolProp*/

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

@@ -199,6 +199,66 @@ namespace SaturationSolvers
{
long double T,p;
};
struct newton_raphson_twophase_options{
enum imposed_variable_options {NO_VARIABLE_IMPOSED = 0, P_IMPOSED, T_IMPOSED};
int Nstep_max;
std::size_t Nsteps;
long double beta, omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap;
imposed_variable_options imposed_variable;
std::vector<long double> x, y, z;
newton_raphson_twophase_options(){ Nstep_max = 30; Nsteps = 0; beta = -1; omega =1;} // Defaults
};
/** \brief A class to do newton raphson solver for VLE for p,q or T,q
*
* A class is used rather than a function so that it is easier to store iteration histories, additional output values, etc.
*
* This class only handles bubble and dew lines since the independent variables are N-1 of the mole fractions in the incipient phase along with one of T, p, or rho
*/
class newton_raphson_twophase
{
public:
newton_raphson_twophase_options::imposed_variable_options imposed_variable;
long double error_rms, rhomolar_liq, rhomolar_vap, T, p, min_rel_change, beta;
unsigned int N;
bool logging;
int Nsteps;
STLMatrix J;
HelmholtzEOSMixtureBackend *HEOS;
std::vector<long double> K, x, y, z, r, negative_r, err_rel;
std::vector<SuccessiveSubstitutionStep> step_logger;
newton_raphson_twophase(){};
void resize(unsigned int N);
// Reset the state of all the internal variables
void pre_call()
{
K.clear(); x.clear(); y.clear(); step_logger.clear(); error_rms = 1e99; Nsteps = 0;
rhomolar_liq = _HUGE; rhomolar_vap = _HUGE; T = _HUGE; p = _HUGE;
};
/** Call the Newton-Raphson VLE Solver
*
* This solver must be passed reasonable guess values for the mole fractions,
* densities, etc. You may want to take a few steps of successive substitution
* before you start with Newton Raphson.
*
* @param HEOS HelmholtzEOSMixtureBackend instance
* @param IO The input/output data structure
*/
void call(HelmholtzEOSMixtureBackend &HEOS, newton_raphson_twophase_options &IO);
/*! Build the arrays for the Newton-Raphson solve
This method builds the Jacobian matrix, the sensitivity matrix, etc.
*
*/
void build_arrays();
};
struct newton_raphson_saturation_options{
enum imposed_variable_options {NO_VARIABLE_IMPOSED = 0, P_IMPOSED, RHOV_IMPOSED, T_IMPOSED};
@@ -208,7 +268,7 @@ namespace SaturationSolvers
long double omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap;
imposed_variable_options imposed_variable;
std::vector<long double> x, y;
newton_raphson_saturation_options(){ Nsteps = 0;} // Defaults
newton_raphson_saturation_options(){ Nstep_max = 30; Nsteps = 0;} // Defaults
};
/** \brief A class to do newton raphson solver for VLE given guess values for vapor-liquid equilibria. This class will then be included in the Mixture class
@@ -226,10 +286,10 @@ namespace SaturationSolvers
bool logging;
bool bubble_point;
int Nsteps;
long double dTsat_dPsat, dPsat_dTsat;
STLMatrix J;
HelmholtzEOSMixtureBackend *HEOS;
std::vector<long double> K, x, y, phi_ij_liq, phi_ij_vap, dlnphi_drho_liq, dlnphi_drho_vap, r, negative_r, dXdS, neg_dFdS, err_rel;
long double dTsat_dPsat, dPsat_dTsat;
std::vector<long double> K, x, y, r, negative_r, err_rel;
std::vector<SuccessiveSubstitutionStep> step_logger;
newton_raphson_saturation(){};
@@ -239,8 +299,7 @@ namespace SaturationSolvers
// Reset the state of all the internal variables
void pre_call()
{
K.clear(); x.clear(); y.clear(); phi_ij_liq.clear();
phi_ij_vap.clear(); dlnphi_drho_liq.clear(), dlnphi_drho_vap.clear(),
K.clear(); x.clear(); y.clear();
step_logger.clear(); error_rms = 1e99; Nsteps = 0;
rhomolar_liq = _HUGE; rhomolar_vap = _HUGE; T = _HUGE; p = _HUGE;
};