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

Newton-raphson works for mixtures

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-30 19:19:10 +02:00
parent d539f5cde8
commit 5e470617c9
3 changed files with 78 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -260,14 +260,34 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
IO.p = io.p;
IO.Nstep_max = 100;
// Dewpoint
NR.call(HEOS, io.y, io.x, IO);
if (std::abs(static_cast<long double>(HEOS._Q)) < 10*DBL_EPSILON){
IO.bubble_point = true;
}
else{
IO.bubble_point = false;
}
IO.x = io.x;
IO.y = io.y;
// 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();
for (;; IO.p *= 1.01)
{
if (IO.bubble_point){
// Bubblepoint, vapor (y) is incipient phase
NR.call(HEOS, IO.x, IO.y, IO);
}
else{
// Dewpoint, liquid (x) is incipient phase
NR.call(HEOS, IO.y, IO.x, IO);
}
std::cout << IO.p << " " << IO.T << " " << IO.rhomolar_liq << " " << IO.rhomolar_vap << std::endl;
// 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();
}
int rr = 4;
}

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

@@ -1104,13 +1104,24 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke
// Reset all the variables and resize
pre_call();
resize(z.size());
y = z;
x = z_incipient;
this->bubble_point = IO.bubble_point;
rhomolar_liq = IO.rhomolar_liq;
rhomolar_vap = IO.rhomolar_vap;
T = IO.T;
p = IO.p;
if (bubble_point){
// Bubblepoint, vapor (y) is the incipient phase
x = z;
y = z_incipient;
}
else{
// Dewpoint, liquid (x) is the incipient phase
y = z;
x = z_incipient;
}
// Hold a pointer to the backend
this->HEOS = &HEOS;
@@ -1123,18 +1134,27 @@ 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::cout << "-r: " << vec_to_string(negative_r, "%0.12Lg") << std::endl;
std::cout << "J: " << vec_to_string(J, "%0.12Lg") << std::endl;
//std::cout << "-r: " << vec_to_string(negative_r, "%0.12Lg") << std::endl;
//std::cout << "J: " << vec_to_string(J, "%0.12Lg") << std::endl;
std::vector<long double> v = linsolve(J, negative_r);
max_rel_change = max_abs_value(v);
for (unsigned int i = 0; i < N-1; ++i){
x[i] += v[i];
}
//std::cout << "v: " << vec_to_string(v, "%0.12Lg") << std::endl;
T += v[N-1];
x[N-1] = 1 - std::accumulate(x.begin(), x.end()-1, 0.0);
std::cout << vec_to_string(x, "%0.12Lg") << std::endl;
if (bubble_point){
for (unsigned int i = 0; i < N-1; ++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){
x[i] += v[i];
}
x[N-1] = 1 - std::accumulate(x.begin(), x.end()-1, 0.0);
}
//std::cout << format("\t%Lg ", this->error_rms) << vec_to_string(v, "%0.10Lg") << " " << vec_to_string(x, "%0.10Lg") << std::endl;
iter++;
}
@@ -1143,6 +1163,9 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke
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;
}
void SaturationSolvers::newton_raphson_saturation::build_arrays()
@@ -1150,7 +1173,6 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
// References to the classes for concision
HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
bool bubble_point = false;
// Step 0:
// -------
// Set mole fractions for the incipient phase
@@ -1164,8 +1186,10 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
}
// Update the liquid and vapor classes
rSatL.update_TP_guessrho(T, p, rhomolar_liq);
rSatV.update_TP_guessrho(T, p, rhomolar_vap);
rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar();
rSatV.update_TP_guessrho(T, p, rhomolar_vap); rhomolar_vap = rSatV.rhomolar();
//std::cout << format("\t\t%Lg %Lg %Lg %s\n", T, rhomolar_liq, rhomolar_vap, vec_to_string(y,"%0.12Lg").c_str());
// For diagnostic purposes calculate the pressures (no derivatives are evaluated)
long double p_liq = rSatL.p();
@@ -1186,9 +1210,19 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
r[i] = ln_f_liq - ln_f_vap;
if (bubble_point){
throw NotImplementedError();
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
if (i == N-1){
J[N-1][j] = -1*(-1/y[N-1] + MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatV,N-1,j, xN_flag));
}
else if (i != j){
J[i][j] = -1*MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatV,i,j, xN_flag);
}
else{
J[i][i] = -1*(1/y[i] + MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatV,i,i, xN_flag));
}
}
}
else{ // its a dewpoint calculation
else{
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
if (i == N-1){
J[N-1][j] = (-1/x[N-1] + MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatL,N-1,j, xN_flag));
@@ -1200,8 +1234,8 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays()
J[i][i] = (1/x[i] + MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatL,i,i, xN_flag));
}
}
J[i][N-1] = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatV, i, xN_flag);
}
J[i][N-1] = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatV, i, xN_flag);
}
// Flip all the signs of the entries in the residual vector since we are solving Jv = -r, not Jv=r

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

@@ -196,6 +196,7 @@ namespace SaturationSolvers
struct newton_raphson_saturation_options{
enum imposed_variable_options {IMPOSED_P, IMPOSED_T};
int Nstep_max;
bool bubble_point;
long double omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T;
imposed_variable_options imposed_variable;
std::vector<long double> x, y;
@@ -214,6 +215,7 @@ namespace SaturationSolvers
long double error_rms, rhomolar_liq, rhomolar_vap, T, p, max_rel_change;
unsigned int N;
bool logging;
bool bubble_point;
int Nsteps;
STLMatrix J;
HelmholtzEOSMixtureBackend *HEOS;