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

Added necessary code to get sat_p working for dewpoint calculations, phase envelopes build for R32/R125 - was getting bad density roots

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-09-11 02:15:19 +02:00
parent 91bedf919a
commit 20cae192b0
7 changed files with 222 additions and 65 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h
View File

@@ -64,6 +64,9 @@ public:
std::vector<long double> &get_K(){return K;};
std::vector<long double> &get_lnK(){return lnK;};
std::vector<long double> calc_mole_fractions_liquid(void){return SatL->get_mole_fractions();};
std::vector<long double> calc_mole_fractions_vapor(void){return SatV->get_mole_fractions();};
const CoolProp::PhaseEnvelopeData &calc_phase_envelope_data(){return PhaseEnvelope;};
void resize(unsigned int N);

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

@@ -96,6 +96,14 @@ long double MixtureDerivatives::d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS,
}
}
long double MixtureDerivatives::dln_fugacity_i_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
return dln_fugacity_coefficient_dT__constp_n(HEOS, i, xN_flag);
}
long double MixtureDerivatives::dln_fugacity_i_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
return dln_fugacity_coefficient_dp__constT_n(HEOS, i, xN_flag) + 1/HEOS.p();
}
long double MixtureDerivatives::fugacity_i(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return HEOS.mole_fractions[i]*HEOS.rhomolar()*HEOS.gas_constant()*HEOS.T()*exp( dnalphar_dni__constT_V_nj(HEOS, i, xN_flag));
}
@@ -148,6 +156,20 @@ long double MixtureDerivatives::dln_fugacity_i_ddelta__consttau_x(HelmholtzEOSMi
{
return 1 + HEOS.delta()*HEOS.dalphar_dDelta() + HEOS.delta()*d_ndalphardni_dDelta(HEOS, i, xN_flag);
}
long double MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
// This is a term to which some more might be added depending on i and j
long double val = dln_fugacity_coefficient_dxj__constT_p_xi(HEOS, i, j, xN_flag);
const std::vector<long double> &x = HEOS.get_const_mole_fractions();
std::size_t N = x.size();
if (i == N-1){
val += -1/x[N-1];
}
else if (i == j){
val += 1/x[j];
}
return val;
}
long double MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
long double rhor = HEOS.Reducing->rhormolar(HEOS.get_const_mole_fractions());
@@ -363,6 +385,8 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
ss00 << "Mixture with " << Ncomp << " components";
SECTION(ss00.str(),"")
{
x_N_dependency_flag xN_flag = XN_INDEPENDENT;
std::vector<std::string> names;
std::vector<long double> z;
shared_ptr<HelmholtzEOSMixtureBackend> HEOS, HEOS_plusT_constrho, HEOS_plusrho_constT, HEOS_minusT_constrho, HEOS_minusrho_constT;
@@ -414,7 +438,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
rHEOS.update(DmolarT_INPUTS, rho1, T1);
x_N_dependency_flag xN_flag = XN_DEPENDENT;
// These ones only require the i index
for (std::size_t i = 0; i< z.size();++i)
{
@@ -423,7 +446,9 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
HEOS_pluszi->specify_phase(iphase_gas);
HelmholtzEOSMixtureBackend &rHEOS_pluszi = *(HEOS_pluszi.get());
std::vector<long double> zp = z; /// Copy base composition
zp[i] += dz; zp[z.size()-1] -= dz;
zp[i] += dz;
if (xN_flag = XN_DEPENDENT)
zp[z.size()-1] -= dz;
rHEOS_pluszi.set_mole_fractions(zp);
rHEOS_pluszi.update(DmolarT_INPUTS, rho1, T1);
@@ -431,7 +456,9 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
HEOS_minuszi->specify_phase(iphase_gas);
HelmholtzEOSMixtureBackend &rHEOS_minuszi = *(HEOS_minuszi.get());
std::vector<long double> zm = z; /// Copy base composition
zm[i] -= dz; zm[z.size()-1] += dz;
zm[i] -= dz;
if (xN_flag = XN_DEPENDENT)
zm[z.size()-1] += dz;
rHEOS_minuszi.set_mole_fractions(zm);
rHEOS_minuszi.update(DmolarT_INPUTS, rho1, T1);
@@ -448,6 +475,19 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
CAPTURE(analytic);
CHECK(err < 1e-8);
}
std::ostringstream ss0a;
ss0a << "dln_fugacity_i_dp__constT_n, i=" << i;
SECTION(ss0a.str(),"")
{
double analytic = MixtureDerivatives::dln_fugacity_i_dp__constT_n(rHEOS, i, xN_flag);
double v1 = log(MixtureDerivatives::fugacity_i(rHEOS_plusrho_constT, i, xN_flag)), p1 = rHEOS_plusrho_constT.p();
double v2 = log(MixtureDerivatives::fugacity_i(rHEOS_minusrho_constT, i, xN_flag)), p2 = rHEOS_minusrho_constT.p();
double numeric = (v1 - v2)/(p1 - p2);
double err = std::abs((numeric-analytic)/analytic);
CAPTURE(numeric);
CAPTURE(analytic);
CHECK(err < 1e-8);
}
std::ostringstream ss0b;
ss0b << "dln_fugacity_i_drho__constT_n, i=" << i;
SECTION(ss0b.str(), "")
@@ -491,6 +531,25 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
}
std::ostringstream ss1a;
ss1a << "dln_fugacity_i_dT__constp_n, i=" << i;
SECTION(ss1a.str(), "")
{
double T1 = 300, dT = 1e-3;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, 101325, T1);
double analytic = MixtureDerivatives::dln_fugacity_i_dT__constp_n(rHEOS, i, xN_flag);
rHEOS.update(PT_INPUTS, 101325, T1 + dT);
double v1 = log(MixtureDerivatives::fugacity_i(rHEOS, i, xN_flag));
rHEOS.update(PT_INPUTS, 101325, T1 - dT);
double v2 = log(MixtureDerivatives::fugacity_i(rHEOS, i, xN_flag));
double numeric = (v1 - v2)/(2*dT);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
}
std::ostringstream ss3;
ss3 << "d_ndalphardni_dDelta, i=" << i;

4
src/Backends/Helmholtz/MixtureDerivatives.h
View File

@@ -124,9 +124,13 @@ class MixtureDerivatives{
*/
static long double dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static long double dln_fugacity_i_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static long double dln_fugacity_i_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static long double dln_fugacity_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static long double dln_fugacity_i_dtau__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static long double dln_fugacity_i_ddelta__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static long double dln_fugacity_dxj__constT_rho_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions
*
* The derivative term

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

@@ -15,7 +15,7 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS)
// Set some imput options
SaturationSolvers::mixture_VLE_IO io;
io.sstype = SaturationSolvers::imposed_p;
io.Nstep_max = 2;
io.Nstep_max = 20;
// Set the pressure to a low pressure
HEOS._p = 100000;
@@ -44,6 +44,22 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS)
IO.p = io.p;
IO.Nstep_max = 30;
/*
IO.p = 1e5;
IO.rhomolar_liq = 17257.17130;
IO.rhomolar_vap = 56.80022884;
IO.T = 219.5200523;
IO.x[0] = 0.6689704673;
IO.x[1] = 0.3310295327;
*/
IO.rhomolar_liq *= 1.2;
IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::P_IMPOSED;
NR.call(HEOS, IO.y, IO.x, IO);
// Switch to density imposed
IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::RHOV_IMPOSED;
bool dont_extrapolate = false;
PhaseEnvelopeData &env = HEOS.PhaseEnvelope;
@@ -63,7 +79,6 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS)
}
if (iter - iter0 > 0){ IO.rhomolar_vap *= factor;}
long double x = IO.rhomolar_vap;
if (dont_extrapolate)
{
// Reset the step to a reasonably small size
@@ -71,29 +86,29 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS)
}
else if (iter - iter0 == 2)
{
IO.T = LinearInterp(env.rhomolar_vap, env.T, iter-2, iter-1, x);
IO.rhomolar_liq = LinearInterp(env.rhomolar_vap, env.rhomolar_liq, iter-2, iter-1, x);
IO.T = LinearInterp(env.rhomolar_vap, env.T, iter-2, iter-1, IO.rhomolar_vap);
IO.rhomolar_liq = LinearInterp(env.rhomolar_vap, env.rhomolar_liq, iter-2, iter-1, IO.rhomolar_vap);
for (std::size_t i = 0; i < IO.x.size()-1; ++i) // First N-1 elements
{
IO.x[i] = LinearInterp(env.rhomolar_vap, env.x[i], iter-2, iter-1, x);
IO.x[i] = LinearInterp(env.rhomolar_vap, env.x[i], iter-2, iter-1, IO.rhomolar_vap);
}
}
else if (iter - iter0 == 3)
{
IO.T = QuadInterp(env.rhomolar_vap, env.T, iter-3, iter-2, iter-1, x);
IO.rhomolar_liq = QuadInterp(env.rhomolar_vap, env.rhomolar_liq, iter-3, iter-2, iter-1, x);
IO.T = QuadInterp(env.rhomolar_vap, env.T, iter-3, iter-2, iter-1, IO.rhomolar_vap);
IO.rhomolar_liq = QuadInterp(env.rhomolar_vap, env.rhomolar_liq, iter-3, iter-2, iter-1, IO.rhomolar_vap);
for (std::size_t i = 0; i < IO.x.size()-1; ++i) // First N-1 elements
{
IO.x[i] = QuadInterp(env.rhomolar_vap, env.x[i], iter-3, iter-2, iter-1, x);
IO.x[i] = QuadInterp(env.rhomolar_vap, env.x[i], iter-3, iter-2, iter-1, IO.rhomolar_vap);
}
}
else if (iter - iter0 > 3)
{
IO.T = CubicInterp(env.rhomolar_vap, env.T, iter-4, iter-3, iter-2, iter-1, x);
IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iter-4, iter-3, iter-2, iter-1, x);
IO.T = CubicInterp(env.rhomolar_vap, env.T, iter-4, iter-3, iter-2, iter-1, IO.rhomolar_vap);
IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, iter-4, iter-3, iter-2, iter-1, IO.rhomolar_vap);
// Check if there is a large deviation from linear interpolation - this suggests a step size that is so large that a minima or maxima of the interpolation function is crossed
long double T_linear = LinearInterp(env.rhomolar_vap, env.T, iter-2, iter-1, x);
long double T_linear = LinearInterp(env.rhomolar_vap, env.T, iter-2, iter-1, IO.rhomolar_vap);
if (std::abs((T_linear-IO.T)/IO.T) > 0.1){
// Try again, but with a smaller step
IO.rhomolar_vap /= factor;
@@ -103,7 +118,7 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS)
}
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], iter-4, iter-3, iter-2, iter-1, x);
IO.x[i] = CubicInterp(env.rhomolar_vap, env.x[i], iter-4, iter-3, iter-2, iter-1, IO.rhomolar_vap);
if (IO.x[i] < 0 || IO.x[i] > 1){
// Try again, but with a smaller step
IO.rhomolar_vap /= factor;
@@ -123,13 +138,18 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS)
// Dewpoint calculation, liquid (x) is incipient phase
try{
NR.call(HEOS, IO.y, IO.x, IO);
if (!ValidNumber(IO.rhomolar_liq) || !ValidNumber(IO.p) || !ValidNumber(IO.T)){
throw ValueError("Invalid number");
}
}
catch(std::exception &e){
std::cout << e.what() << std::endl;
std::cout << IO.T << " " << IO.p << std::endl;
// Try again, but with a smaller step
IO.rhomolar_vap /= factor;
factor = 1 + (factor-1)/2;
failure_count++;
continue;
}

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

@@ -722,7 +722,7 @@ void SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend &HEOS
long double rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(T, p, iphase_liquid); // [mol/m^3]
long double rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(T, p, iphase_gas); // [mol/m^3]
HEOS.SatL->update_TP_guessrho(T, p, rhomolar_liq*1.3);
HEOS.SatL->update_TP_guessrho(T, p, rhomolar_liq*1.5);
HEOS.SatV->update_TP_guessrho(T, p, rhomolar_vap);
do
@@ -796,9 +796,24 @@ void SaturationSolvers::newton_raphson_saturation::resize(unsigned int N)
x.resize(N);
y.resize(N);
r.resize(N+1);
negative_r.resize(N+1);
J.resize(N+1, std::vector<long double>(N+1, 0));
if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){
r.resize(N+1);
negative_r.resize(N+1);
J.resize(N+1);
for (std::size_t i = 0; i < N+1; ++i){
J[i].resize(N+1);
}
err_rel.resize(N+1);
}
else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED){
r.resize(N);
negative_r.resize(N);
J.resize(N, std::vector<long double>(N, 0));
err_rel.resize(N);
}
else{
throw ValueError();
}
}
void SaturationSolvers::newton_raphson_saturation::check_Jacobian()
{
@@ -893,16 +908,18 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke
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();
resize(z.size());
this->bubble_point = IO.bubble_point;
rhomolar_liq = IO.rhomolar_liq;
rhomolar_vap = IO.rhomolar_vap;
T = IO.T;
p = IO.p;
imposed_variable = IO.imposed_variable;
resize(z.size());
if (bubble_point){
// Bubblepoint, vapor (y) is the incipient phase
@@ -928,9 +945,6 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke
// Solve for the step; v is the step with the contents
// [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)]
std::vector<long double> v = linsolve(J, negative_r);
max_rel_change = max_abs_value(v);
min_abs_change = min_abs_value(v);
if (bubble_point){
for (unsigned int i = 0; i < N-1; ++i){
@@ -940,21 +954,26 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke
}
else{
for (unsigned int i = 0; i < N-1; ++i){
err_rel[i] = v[i]/y[i];
x[i] += v[i];
}
x[N-1] = 1 - std::accumulate(x.begin(), x.end()-1, 0.0);
}
T += v[N-1];
rhomolar_liq += v[N];
//std::cout << format("\t%Lg ", this->error_rms) << T << " " << rhomolar_liq << " " << rhomolar_vap << " " << vec_to_string(v, "%0.10Lg") << " " << vec_to_string(x, "%0.10Lg") << std::endl;
T += v[N-1]; err_rel[N-1] = v[N-1]/T;
if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){
rhomolar_liq += v[N]; err_rel[N] = v[N]/rhomolar_liq;
}
//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-7 && max_rel_change > 1000*LDBL_EPSILON && min_abs_change > 100*DBL_EPSILON && iter < IO.Nstep_max);
while(this->error_rms > 1e-9 && min_rel_change > 10*DBL_EPSILON && iter < IO.Nstep_max);
IO.Nsteps = iter;
IO.p = p;
IO.x = x; // Mole fractions in liquid
@@ -983,54 +1002,93 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
else{
// Liquid is incipient phase, set its composition
rSatL.set_mole_fractions(x);
rSatV.set_mole_fractions(y);
}
rSatL.update(DmolarT_INPUTS, rhomolar_liq, T);
rSatV.update(DmolarT_INPUTS, rhomolar_vap, T);
//std::cout << format("\t\t%Lg %Lg %Lg %s\n", T, rhomolar_liq, rhomolar_vap, vec_to_string(x,"%0.12Lg").c_str());
if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){
rSatL.update(DmolarT_INPUTS, rhomolar_liq, T);
rSatV.update(DmolarT_INPUTS, rhomolar_vap, T);
}
else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED){
rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar();
rSatV.update_TP_guessrho(T, p, rhomolar_vap);
}
else{
throw ValueError();
}
// 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;
// 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 (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED){
x_N_dependency_flag xN_flag = XN_DEPENDENT;
// 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_rho_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_rho_xi(rSatL, i, j, xN_flag);
}
}
J[i][N-1] = MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatV, i, xN_flag);
J[i][N] = MixtureDerivatives::dln_fugacity_i_drho__constT_n(rSatL, i, xN_flag);
}
// ---------------------------------------------------------------
// Derivatives of pL(T,rho',x)-p(T,rho'',y) with respect to inputs
// ---------------------------------------------------------------
r[N] = p_liq - p_vap;
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
J[N][j] = MixtureDerivatives::dpdxj__constT_V_xi(rSatL, j, xN_flag); // p'' not a function of x0
}
// Fixed composition derivatives
J[N][N-1] = rSatL.first_partial_deriv(iP, iT, iDmolar)-rSatV.first_partial_deriv(iP, iT, iDmolar);
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'
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_rho_xi(rSatV, i, j, xN_flag);
// 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_rho_xi(rSatL, 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_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatV, i, xN_flag);
}
J[i][N-1] = MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatV, i, xN_flag);
J[i][N] = MixtureDerivatives::dln_fugacity_i_drho__constT_n(rSatL, i, xN_flag);
}
// ---------------------------------------------------------------
// Derivatives of pL(T,rho',x)-p(T,rho'',y) with respect to inputs
// ---------------------------------------------------------------
r[N] = p_liq - p_vap;
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
J[N][j] = MixtureDerivatives::dpdxj__constT_V_xi(rSatL, j, xN_flag); // p'' not a function of x0
else{
throw ValueError();
}
// Fixed composition derivatives
J[N][N-1] = rSatL.first_partial_deriv(iP, iT, iDmolar)-rSatV.first_partial_deriv(iP, iT, iDmolar);
J[N][N] = rSatL.first_partial_deriv(iP, iDmolar, iT);
// 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

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

@@ -201,7 +201,7 @@ namespace SaturationSolvers
};
struct newton_raphson_saturation_options{
enum imposed_variable_options {IMPOSED_P, IMPOSED_T};
enum imposed_variable_options {NO_VARIABLE_IMPOSED = 0, P_IMPOSED, RHOV_IMPOSED};
int Nstep_max;
bool bubble_point;
std::size_t Nsteps;
@@ -219,8 +219,9 @@ namespace SaturationSolvers
*/
class newton_raphson_saturation
{
public:
long double error_rms, rhomolar_liq, rhomolar_vap, T, p, max_rel_change, min_abs_change;
public:
newton_raphson_saturation_options::imposed_variable_options imposed_variable;
long double error_rms, rhomolar_liq, rhomolar_vap, T, p, min_rel_change;
unsigned int N;
bool logging;
bool bubble_point;
@@ -228,7 +229,7 @@ namespace SaturationSolvers
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;
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;
std::vector<SuccessiveSubstitutionStep> step_logger;
newton_raphson_saturation(){};