github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-04-23 03:00:17 -04:00
Code Issues Packages Projects Releases Wiki Activity

First cut at new NR solver for bubble/dew points with mixtures. Sort of works, seems like there is a problem with one of the derivatives

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-26 13:47:34 +02:00
parent 07cc87b8b2
commit e0e0b1e234
5 changed files with 308 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -252,16 +252,16 @@ 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);
//SaturationSolvers::newton_raphson_saturation NR;
//SaturationSolvers::newton_raphson_saturation_options IO;
//IO.rhomolar_liq = io.rhomolar_liq;
//IO.rhomolar_vap = io.rhomolar_vap;
//IO.T = io.T;
//IO.p = io.p;
//IO.Nstep_max = 25;
SaturationSolvers::newton_raphson_saturation NR;
SaturationSolvers::newton_raphson_saturation_options IO;
IO.rhomolar_liq = io.rhomolar_liq;
IO.rhomolar_vap = io.rhomolar_vap;
IO.T = io.T;
IO.p = io.p;
IO.Nstep_max = 25;
// Dewpoint
//NR.call(HEOS, io.y, io.x, IO);
NR.call(HEOS, io.y, io.x, IO);
// Load the outputs
HEOS._phase = iphase_twophase;

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

@@ -19,6 +19,9 @@ long double MixtureDerivatives::d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS,
return 0 + HEOS.Excess.d2alphardxidxj(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i, j);
}
long double MixtureDerivatives::fugacity_i(HelmholtzEOSMixtureBackend &HEOS, std::size_t i){
return HEOS.mole_fractions[i]*HEOS.rhomolar()*HEOS.gas_constant()*HEOS.T()*exp( dnalphar_dni__constT_V_nj(HEOS, i));
}
long double MixtureDerivatives::ln_fugacity_coefficient(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
{
return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i)-log(1+HEOS._delta.pt()*HEOS.dalphar_dDelta());
@@ -68,7 +71,7 @@ long double MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(Helmho
// Gernert 3.115
long double R_u = HEOS.gas_constant();
// partial molar volume is -dpdn/dpdV, so need to flip the sign here
return d2nalphar_dxi_dnj__constT_V(HEOS, i, j) - partial_molar_volume(HEOS, i)/(R_u*HEOS._T)*dpdxj__constT_V_xi(HEOS, j);
return d2nalphar_dxj_dni__constT_V(HEOS, j, i) - partial_molar_volume(HEOS, i)/(R_u*HEOS._T)*dpdxj__constT_V_xi(HEOS, j);
}
long double MixtureDerivatives::dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j)
{

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

@@ -76,6 +76,15 @@ class MixtureDerivatives{
*/
static long double partial_molar_volume(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
/** \brief Fugacity of the i-th component
*
* Given by the equation
* \f[
* f_i(\delta, \tau, \bar x) = x_i\rho R T \exp\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_{j \neq i}}
* \f]
*/
static long double fugacity_i(HelmholtzEOSMixtureBackend &HEOS, std::size_t i);
/** \brief Natural logarithm of the fugacity coefficient
*
* @param HEOS The HelmholtzEOSMixtureBackend to be used
@@ -172,11 +181,11 @@ class MixtureDerivatives{
*
* The derivative term
* \f[
* \left(\frac{\partial^2n\alpha^r}{\partial x_i\partial n_j} \right)_{T,V} = \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_{j\neq i}}\right)\right)_{T,V,x_{k\neq j}} +\left(\frac{\partial \alpha^r}{\partial x_j}\right)_{T,V,x_{k\neq j}}
* \left(\frac{\partial^2n\alpha^r}{\partial x_j\partial n_i} \right)_{T,V} = \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_{j\neq i}}\right)\right)_{T,V,x_{k\neq j}} +\left(\frac{\partial \alpha^r}{\partial x_j}\right)_{T,V,x_{k\neq j}}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double d2nalphar_dxi_dnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j){ return MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HEOS, i, j) + MixtureDerivatives::dalphar_dxj__constT_V_xi(HEOS, j);};
static long double d2nalphar_dxj_dni__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, std::size_t i){ return MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HEOS, i, j) + MixtureDerivatives::dalphar_dxj__constT_V_xi(HEOS, j);};
/// Gernert Equation 3.119
/// Catch test provided

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

@@ -790,6 +790,8 @@ void SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend &HEOS
options.T = HEOS.SatL->T();
options.rhomolar_liq = HEOS.SatL->rhomolar();
options.rhomolar_vap = HEOS.SatV->rhomolar();
options.x = x;
options.y = y;
}
void SaturationSolvers::newton_raphson_VLE_GV::resize(unsigned int N)
@@ -879,7 +881,6 @@ void SaturationSolvers::newton_raphson_VLE_GV::check_Jacobian(HelmholtzEOSMixtur
}
std::cout << vec_to_string(get_col(Jnum, N+1),"%12.11f") << std::endl;
std::cout << vec_to_string(get_col(J,N+1),"%12.11f") << std::endl;
}
void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &K, mixture_VLE_IO &IO)
{
@@ -930,10 +931,9 @@ void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend &
while(this->error_rms > 1e-8 && max_rel_change > 1000*LDBL_EPSILON && iter < IO.Nstep_max);
Nsteps = iter;
IO.p = p;
IO.x = &x; // Mole fractions in liquid
IO.y = &y; // Mole fractions in vapor
IO.x = x; // Mole fractions in liquid
IO.y = y; // Mole fractions in vapor
}
void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureBackend &HEOS, long double beta, long double T, long double rhomolar_liq, const long double rhomolar_vap, const std::vector<long double> &z, std::vector<long double> &K)
{
// Step 0:
@@ -1032,6 +1032,221 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB
error_rms = sqrt(error_rms); // Square-root (The R in RMS)
}
void SaturationSolvers::newton_raphson_saturation::resize(unsigned int N)
{
this->N = N;
x.resize(N);
y.resize(N);
r.resize(N);
negative_r.resize(N);
J.resize(N, std::vector<long double>(N, 0));
}
void SaturationSolvers::newton_raphson_saturation::check_Jacobian()
{
// References to the classes for concision
HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
// Build the Jacobian and residual vectors
build_arrays();
// Make copies of the base
long double T0 = T;
std::vector<long double> r0 = r, x0 = x;
STLMatrix J0 = J;
// Derivatives with respect to T
double dT = 1e-3, T1 = T+dT, T2 = T-dT;
T = T1;
rSatL.update_TP_guessrho(T,p,rhomolar_liq); rhomolar_liq = rSatL.rhomolar();
rSatV.update_TP_guessrho(T,p,rhomolar_vap); rhomolar_vap = rSatV.rhomolar();
build_arrays();
std::vector<long double> r1 = r;
T = T2;
rSatL.update_TP_guessrho(T,p,rhomolar_liq); rhomolar_liq = rSatL.rhomolar();
rSatV.update_TP_guessrho(T,p,rhomolar_vap); rhomolar_vap = rSatV.rhomolar();
build_arrays();
std::vector<long double> r2 = r;
std::vector<long double> diffn(N, _HUGE);
for (std::size_t i = 0; i < N; ++i){
diffn[i] = (r1[i]-r2[i])/(2*dT);
}
std::cout << "numerical: " << vec_to_string(diffn, "%0.11Lg") << std::endl;
std::cout << "analytic: " << vec_to_string(get_col(J0, N-1), "%0.11Lg") << std::endl;
// Derivatives with respect to x0
double dx = 1e-5, T = T0;
x = x0; x[0] += dx;// x[1] -= dx;
rSatL.set_mole_fractions(x);
rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar();
build_arrays(); r1 = r;
x = x0; x[0] -= dx;// x[1] += dx;
rSatL.set_mole_fractions(x);
rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar();
build_arrays(); r2 = r;
for (std::size_t i = 0; i < N; ++i){
diffn[i] = (r1[i]-r2[i])/(2*dx);
}
std::cout << "numerical: " << vec_to_string(diffn, "%0.11Lg") << std::endl;
std::cout << "analytic: " << vec_to_string(get_col(J0, 0), "%0.11Lg") << std::endl;
int rrr = 0;
// Jnum.resize(N, std::vector<long double>(N, 0));
//
// for (std::size_t j = 0; j < N; ++j)
// {
// std::vector<long double> KK = K;
// KK[j] += 1e-6;
// build_arrays(HEOS,beta0,T0,rhomolar_liq0,rhomolar_vap0,z,KK);
// std::vector<long double> r1 = r;
// for (std::size_t i = 0; i < N+2; ++i)
// {
// Jnum[i][j] = (r1[i]-r0[i])/(log(KK[j])-log(K[j]));
// }
// std::cout << vec_to_string(get_col(Jnum,j),"%12.11f") << std::endl;
// std::cout << vec_to_string(get_col(J,j),"%12.11f") << std::endl;
// }
//
// build_arrays(HEOS,beta0,T0+1e-6,rhomolar_liq0,rhomolar_vap0,z,K);
// std::vector<long double> r1 = r, JN = r;
// for (std::size_t i = 0; i < r.size(); ++i)
// {
// Jnum[i][N] = (r1[i]-r0[i])/(log(T0+1e-6)-log(T0));
// }
// std::cout << vec_to_string(get_col(Jnum,N),"%12.11f") << std::endl;
// std::cout << vec_to_string(get_col(J,N),"%12.11f") << std::endl;
//
// // Build the Jacobian and residual vectors for given inputs of K_i,T,p
// build_arrays(HEOS,beta0,T0,rhomolar_liq0+1e-3,rhomolar_vap0,z,K);
// std::vector<long double> r2 = r, JNp1 = r;
// for (std::size_t i = 0; i < r.size(); ++i)
// {
// Jnum[i][N+1] = (r2[i]-r0[i])/(log(rhomolar_liq0+1e-3)-log(rhomolar_liq0));
// }
// std::cout << vec_to_string(get_col(Jnum, N+1),"%12.11f") << std::endl;
// std::cout << vec_to_string(get_col(J,N+1),"%12.11f") << std::endl;
}
void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &z_incipient, newton_raphson_saturation_options &IO)
{
int iter = 0;
// Reset all the variables and resize
pre_call();
resize(z.size());
y = z;
x = z_incipient;
rhomolar_liq = IO.rhomolar_liq;
rhomolar_vap = IO.rhomolar_vap;
T = IO.T;
p = IO.p;
// Hold a pointer to the backend
this->HEOS = &HEOS;
check_Jacobian();
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)]
std::cout << vec_to_string(negative_r, "%0.12Lg") << std::endl;
std::cout << 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];
}
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;
iter++;
}
while(this->error_rms > 1e-12 && max_rel_change > 1000*LDBL_EPSILON && iter < IO.Nstep_max);
Nsteps = iter;
IO.p = p;
IO.x = x; // Mole fractions in liquid
IO.y = y; // Mole fractions in vapor
}
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
if (bubble_point){
// Vapor is incipient phase, set its composition
rSatV.set_mole_fractions(y);
}
else{
// Liquid is incipient phase, set its composition
rSatL.set_mole_fractions(x);
}
// Update the liquid and vapor classes
rSatL.update_TP_guessrho(T, p, rhomolar_liq);
rSatV.update_TP_guessrho(T, p, rhomolar_vap);
// 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
// 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));
long double ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i));
r[i] = ln_f_liq - ln_f_vap;
if (bubble_point){
throw NotImplementedError();
}
else{
for (std::size_t j = 0; j < N-1; ++j){ // j from 0 to N-2
if (i == N-1){
J[i][j] = (-1/x[i] + MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatL,i,j));
}
else if (i != j){
J[i][j] = MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatL,i,j);
}
else{
J[i][i] = (1/x[i] + MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rSatL,i,i));
}
}
J[i][N-1] = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatL, i) - MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatV, i);
}
}
// 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 < N; ++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)
}
void PhaseEnvelope::PhaseEnvelope_GV::build(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &K, SaturationSolvers::mixture_VLE_IO &IO)
{
// Use the residual function based on ln(K_i), ln(T) and ln(rho') as independent variables. rho'' is specified
@@ -1083,7 +1298,7 @@ void PhaseEnvelope::PhaseEnvelope_GV::build(HelmholtzEOSMixtureBackend &HEOS, co
NRVLE.check_Jacobian(HEOS,z,K,IO_NRVLE);
}*/
NRVLE.call(HEOS,z,K,IO_NRVLE);
dew.store_variables(IO_NRVLE.T,IO_NRVLE.p,IO_NRVLE.rhomolar_liq,IO_NRVLE.rhomolar_vap,K,*(IO_NRVLE.x),*(IO_NRVLE.y));
dew.store_variables(IO_NRVLE.T,IO_NRVLE.p,IO_NRVLE.rhomolar_liq,IO_NRVLE.rhomolar_vap,K,IO_NRVLE.x,IO_NRVLE.y);
iter ++;
std::cout << format("%g %g %g %g %g %d %g\n",IO_NRVLE.p,IO_NRVLE.rhomolar_liq,IO_NRVLE.rhomolar_vap,IO_NRVLE.T,K[0],NRVLE.Nsteps,factor);
if (iter < 5){continue;}

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

@@ -33,7 +33,7 @@ namespace SaturationSolvers
{
int sstype, Nstep_max;
long double rhomolar_liq, rhomolar_vap, p, T, beta;
std::vector<long double> *x, *y, *K;
std::vector<long double> x, y, K;
};
/*! Returns the natural logarithm of K for component i using the method from Wilson as in
@@ -193,6 +193,70 @@ namespace SaturationSolvers
long double T,p;
};
struct newton_raphson_saturation_options{
enum imposed_variable_options {IMPOSED_P, IMPOSED_T};
int Nstep_max;
long double omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T;
imposed_variable_options imposed_variable;
std::vector<long double> x, y;
newton_raphson_saturation_options(){ } // 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
*
* 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_saturation
{
public:
long double error_rms, rhomolar_liq, rhomolar_vap, T, p, max_rel_change;
unsigned int N;
bool logging;
int Nsteps;
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<SuccessiveSubstitutionStep> step_logger;
newton_raphson_saturation(){};
void resize(unsigned int N);
// 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(),
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 z Bulk mole fractions [-]
@param z_incipient Initial guesses for the mole fractions of the incipient phase [-]
*/
void call(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &z_incipient, newton_raphson_saturation_options &IO);
/*! Build the arrays for the Newton-Raphson solve
This method builds the Jacobian matrix, the sensitivity matrix, etc.
*
*/
void build_arrays();
/** Check the derivatives in the Jacobian using numerical derivatives.
*/
void check_Jacobian();
};
/*!
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