From 5e470617c92ee0001cdbbfc295cc454a42a3e022 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Sat, 30 Aug 2014 19:19:10 +0200 Subject: [PATCH] Newton-raphson works for mixtures Signed-off-by: Ian Bell --- src/Backends/Helmholtz/FlashRoutines.cpp | 34 ++++++++++--- src/Backends/Helmholtz/VLERoutines.cpp | 64 ++++++++++++++++++------ src/Backends/Helmholtz/VLERoutines.h | 2 + 3 files changed, 78 insertions(+), 22 deletions(-) diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index dcda1e5f..0a606a6c 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -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(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; } diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index c67b5097..1384054b 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -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 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 diff --git a/src/Backends/Helmholtz/VLERoutines.h b/src/Backends/Helmholtz/VLERoutines.h index fc12bb99..48934016 100644 --- a/src/Backends/Helmholtz/VLERoutines.h +++ b/src/Backends/Helmholtz/VLERoutines.h @@ -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 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;