CoolProp 6.8.1dev
An open-source fluid property and humid air property database
HelmholtzEOSMixtureBackend.cpp
Go to the documentation of this file.
1/*
2 * AbstractBackend.cpp
3 *
4 * Created on: 20 Dec 2013
5 * Author: jowr
6 */
7
8#include <memory>
9
10#if defined(_MSC_VER)
11# define _CRTDBG_MAP_ALLOC
12# ifndef _CRT_SECURE_NO_WARNINGS
13# define _CRT_SECURE_NO_WARNINGS
14# endif
15# include <crtdbg.h>
16# include <sys/stat.h>
17#else
18# include <sys/stat.h>
19#endif
20
21#include <string>
22//#include "CoolProp.h"
23
24#include "HelmholtzEOSMixtureBackend.h"
25#include "HelmholtzEOSBackend.h"
26#include "Fluids/FluidLibrary.h"
27#include "Solvers.h"
28#include "MatrixMath.h"
29#include "VLERoutines.h"
30#include "FlashRoutines.h"
31#include "TransportRoutines.h"
32#include "MixtureDerivatives.h"
33#include "PhaseEnvelopeRoutines.h"
34#include "ReducingFunctions.h"
35#include "MixtureParameters.h"
36#include "IdealCurves.h"
37#include "MixtureParameters.h"
38#include <stdlib.h>
39
40static int deriv_counter = 0;
41
42namespace CoolProp {
43
44class HEOSGenerator : public AbstractStateGenerator
45{
46 public:
47 AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
48 if (fluid_names.size() == 1) {
49 return new HelmholtzEOSBackend(fluid_names[0]);
50 } else {
51 return new HelmholtzEOSMixtureBackend(fluid_names);
52 }
53 };
54};
55// This static initialization will cause the generator to register
56static CoolProp::GeneratorInitializer<HEOSGenerator> heos_gen(CoolProp::HEOS_BACKEND_FAMILY);
57
58HelmholtzEOSMixtureBackend::HelmholtzEOSMixtureBackend() {
59 imposed_phase_index = iphase_not_imposed;
60 is_pure_or_pseudopure = false;
61 N = 0;
62 _phase = iphase_unknown;
63 // Reset the residual Helmholtz energy class
64 residual_helmholtz.reset(new ResidualHelmholtz());
65}
66HelmholtzEOSMixtureBackend::HelmholtzEOSMixtureBackend(const std::vector<std::string>& component_names, bool generate_SatL_and_SatV) {
67 std::vector<CoolPropFluid> components(component_names.size());
68 for (unsigned int i = 0; i < components.size(); ++i) {
69 components[i] = get_library().get(component_names[i]);
70 }
71
72 // Reset the residual Helmholtz energy class
73 residual_helmholtz.reset(new ResidualHelmholtz());
74
75 // Set the components and associated flags
76 set_components(components, generate_SatL_and_SatV);
77
78 // Set the phase to default unknown value
79 _phase = iphase_unknown;
80}
81HelmholtzEOSMixtureBackend::HelmholtzEOSMixtureBackend(const std::vector<CoolPropFluid>& components, bool generate_SatL_and_SatV) {
82
83 // Reset the residual Helmholtz energy class
84 residual_helmholtz.reset(new ResidualHelmholtz());
85
86 // Set the components and associated flags
87 set_components(components, generate_SatL_and_SatV);
88
89 // Set the phase to default unknown value
90 _phase = iphase_unknown;
91}
92void HelmholtzEOSMixtureBackend::set_components(const std::vector<CoolPropFluid>& components, bool generate_SatL_and_SatV) {
93
94 // Copy the components
95 this->components = components;
96 this->N = components.size();
97
98 is_pure_or_pseudopure = (components.size() == 1);
99 if (is_pure_or_pseudopure) {
100 mole_fractions = std::vector<CoolPropDbl>(1, 1);
101 std::vector<std::vector<double>> ones(1, std::vector<double>(1, 1));
102 Reducing = shared_ptr<ReducingFunction>(new GERG2008ReducingFunction(components, ones, ones, ones, ones));
103 _reducing = calc_reducing_state_nocache(mole_fractions);
104 _gas_constant = calc_gas_constant();
105 } else {
106 // Set the mixture parameters - binary pair reducing functions, departure functions, F_ij, etc.
107 set_mixture_parameters();
108 }
109
110 imposed_phase_index = iphase_not_imposed;
111
112 // Top-level class can hold copies of the base saturation classes,
113 // saturation classes cannot hold copies of the saturation classes
114 if (generate_SatL_and_SatV) {
115 SatL.reset(get_copy(false));
116 SatL->specify_phase(iphase_liquid);
117 linked_states.push_back(SatL);
118 SatL->clear();
119 SatV.reset(get_copy(false));
120 SatV->specify_phase(iphase_gas);
121 SatV->clear();
122 linked_states.push_back(SatV);
123 }
124}
125void HelmholtzEOSMixtureBackend::set_mole_fractions(const std::vector<CoolPropDbl>& mf) {
126 if (mf.size() != N) {
127 throw ValueError(format("size of mole fraction vector [%d] does not equal that of component vector [%d]", mf.size(), N));
128 }
129 // Copy values without reallocating memory
130 this->mole_fractions = mf; // Most effective copy
131 this->resize(N); // No reallocation of this->mole_fractions happens
132 clear_comp_change();
133};
134void HelmholtzEOSMixtureBackend::sync_linked_states(const HelmholtzEOSMixtureBackend* const source) {
135 residual_helmholtz.reset(source->residual_helmholtz->copy_ptr());
136 if (source->Reducing) {
137 Reducing.reset(source->Reducing->copy());
138 }
139 // Recurse into linked states of the class
140 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
141 it->get()->sync_linked_states(source);
142 }
143}
144HelmholtzEOSMixtureBackend* HelmholtzEOSMixtureBackend::get_copy(bool generate_SatL_and_SatV) {
145 // Set up the class with these components
146 HelmholtzEOSMixtureBackend* ptr = new HelmholtzEOSMixtureBackend(components, generate_SatL_and_SatV);
147 // Recursively walk into linked states, setting the departure and reducing terms
148 // to be equal to the parent (this instance)
149 ptr->sync_linked_states(this);
150 return ptr;
151};
152void HelmholtzEOSMixtureBackend::set_mass_fractions(const std::vector<CoolPropDbl>& mass_fractions) {
153 if (mass_fractions.size() != N) {
154 throw ValueError(format("size of mass fraction vector [%d] does not equal that of component vector [%d]", mass_fractions.size(), N));
155 }
156 std::vector<CoolPropDbl> moles;
157 CoolPropDbl sum_moles = 0.0;
158 CoolPropDbl tmp = 0.0;
159 for (unsigned int i = 0; i < components.size(); ++i) {
160 tmp = mass_fractions[i] / components[i].molar_mass();
161 moles.push_back(tmp);
162 sum_moles += tmp;
163 }
164 std::vector<CoolPropDbl> mole_fractions;
165 for (std::vector<CoolPropDbl>::iterator it = moles.begin(); it != moles.end(); ++it) {
166 mole_fractions.push_back(*it / sum_moles);
167 }
168 this->set_mole_fractions(mole_fractions);
169};
170void HelmholtzEOSMixtureBackend::resize(std::size_t N) {
171 this->mole_fractions.resize(N);
172 this->K.resize(N);
173 this->lnK.resize(N);
174 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
175 it->get()->N = N;
176 it->get()->resize(N);
177 }
178}
179void HelmholtzEOSMixtureBackend::recalculate_singlephase_phase() {
180 if (p() > p_critical()) {
181 if (T() > T_critical()) {
182 _phase = iphase_supercritical;
183 } else {
184 _phase = iphase_supercritical_liquid;
185 }
186 } else {
187 if (T() > T_critical()) {
188 _phase = iphase_supercritical_gas;
189 } else {
190 // Liquid or vapor
191 if (rhomolar() > rhomolar_critical()) {
192 _phase = iphase_liquid;
193 } else {
194 _phase = iphase_gas;
195 }
196 }
197 }
198}
199std::string HelmholtzEOSMixtureBackend::fluid_param_string(const std::string& ParamName) {
200 CoolProp::CoolPropFluid cpfluid = get_components()[0];
201 if (!ParamName.compare("name")) {
202 return cpfluid.name;
203 } else if (!ParamName.compare("aliases")) {
204 return strjoin(cpfluid.aliases, get_config_string(LIST_STRING_DELIMITER));
205 } else if (!ParamName.compare("CAS") || !ParamName.compare("CAS_number")) {
206 return cpfluid.CAS;
207 } else if (!ParamName.compare("formula")) {
208 return cpfluid.formula;
209 } else if (!ParamName.compare("ASHRAE34")) {
210 return cpfluid.environment.ASHRAE34;
211 } else if (!ParamName.compare("REFPROPName") || !ParamName.compare("REFPROP_name") || !ParamName.compare("REFPROPname")) {
212 return cpfluid.REFPROPname;
213 } else if (ParamName.find("BibTeX") == 0) // Starts with "BibTeX"
214 {
215 std::vector<std::string> parts = strsplit(ParamName, '-');
216 if (parts.size() != 2) {
217 throw ValueError(format("Unable to parse BibTeX string %s", ParamName.c_str()));
218 }
219 std::string key = parts[1];
220 if (!key.compare("EOS")) {
221 return cpfluid.EOS().BibTeX_EOS;
222 } else if (!key.compare("CP0")) {
223 return cpfluid.EOS().BibTeX_CP0;
224 } else if (!key.compare("VISCOSITY")) {
225 return cpfluid.transport.BibTeX_viscosity;
226 } else if (!key.compare("CONDUCTIVITY")) {
227 return cpfluid.transport.BibTeX_conductivity;
228 } else if (!key.compare("ECS_LENNARD_JONES")) {
229 throw NotImplementedError();
230 } else if (!key.compare("ECS_VISCOSITY_FITS")) {
231 throw NotImplementedError();
232 } else if (!key.compare("ECS_CONDUCTIVITY_FITS")) {
233 throw NotImplementedError();
234 } else if (!key.compare("SURFACE_TENSION")) {
235 return cpfluid.ancillaries.surface_tension.BibTeX;
236 } else if (!key.compare("MELTING_LINE")) {
237 return cpfluid.ancillaries.melting_line.BibTeX;
238 } else {
239 throw CoolProp::KeyError(format("Bad key to get_BibTeXKey [%s]", key.c_str()));
240 }
241 } else if (ParamName.find("pure") == 0) {
242 if (is_pure()) {
243 return "true";
244 } else {
245 return "false";
246 }
247 } else if (ParamName == "INCHI" || ParamName == "InChI" || ParamName == "INCHI_STRING") {
248 return cpfluid.InChI;
249 } else if (ParamName == "INCHI_Key" || ParamName == "InChIKey" || ParamName == "INCHIKEY") {
250 return cpfluid.InChIKey;
251 } else if (ParamName == "2DPNG_URL") {
252 return cpfluid.TwoDPNG_URL;
253 } else if (ParamName == "SMILES" || ParamName == "smiles") {
254 return cpfluid.smiles;
255 } else if (ParamName == "CHEMSPIDER_ID") {
256 return format("%d", cpfluid.ChemSpider_id);
257 } else if (ParamName == "JSON") {
258 return get_fluid_as_JSONstring(cpfluid.CAS);
259 } else {
260 throw ValueError(format("fluid parameter [%s] is invalid", ParamName.c_str()));
261 }
262}
263
264std::optional<EquationOfState::SuperAncillary_t>& HelmholtzEOSMixtureBackend::get_superanc_optional(){
265 if (!is_pure()){
266 throw CoolProp::ValueError("Only available for pure (and not pseudo-pure) fluids");
267 }
268 return components[0].EOS().get_superanc_optional();
269}
270
271double HelmholtzEOSMixtureBackend::get_fluid_parameter_double(const size_t i, const std::string& parameter){
272 if (i >= N) {
273 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N-1));
274 }
275 auto& optsuperanc = get_superanc_optional();
276 if (parameter.find("SUPERANC::") == 0){
277 auto& superanc = optsuperanc.value();
278 if (optsuperanc){
279 std::string key = parameter.substr(10);
280 if (key == "pmax"){
281 return superanc.get_pmax();
282 }
283 else if (key == "pmin"){
284 return superanc.get_pmin();
285 }
286 else if (key == "Tmin"){
287 return superanc.get_Tmin();
288 }
289 else if (key == "Tcrit_num"){
290 return superanc.get_Tcrit_num();
291 }
292 else {
293 throw ValueError(format("Superancillary parameter [%s] is invalid", key.c_str()));
294 }
295 }
296 else{
297 throw ValueError(format("Superancillary not available for this fluid"));
298 }
299 } else {
300 throw ValueError(format("fluid parameter [%s] is invalid", parameter.c_str()));
301 }
302}
303
304
305void HelmholtzEOSMixtureBackend::apply_simple_mixing_rule(std::size_t i, std::size_t j, const std::string& model) {
306 // bound-check indices
307 if (i < 0 || i >= N) {
308 if (j < 0 || j >= N) {
309 throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N-1));
310 } else {
311 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N-1));
312 }
313 } else if (j < 0 || j >= N) {
314 throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, N-1));
315 }
316 if (model == "linear") {
317 double Tc1 = get_fluid_constant(i, iT_critical), Tc2 = get_fluid_constant(j, iT_critical);
318 double gammaT = 0.5 * (Tc1 + Tc2) / sqrt(Tc1 * Tc2);
319 double rhoc1 = get_fluid_constant(i, irhomolar_critical), rhoc2 = get_fluid_constant(j, irhomolar_critical);
320 double gammaV = 4.0 * (1 / rhoc1 + 1 / rhoc2) / pow(pow(rhoc1, -1.0 / 3.0) + pow(rhoc2, -1.0 / 3.0), 3);
321 set_binary_interaction_double(i, j, "betaT", 1.0);
322 set_binary_interaction_double(i, j, "gammaT", gammaT);
323 set_binary_interaction_double(i, j, "betaV", 1.0);
324 set_binary_interaction_double(i, j, "gammaV", gammaV);
325 } else if (model == "Lorentz-Berthelot") {
326 set_binary_interaction_double(i, j, "betaT", 1.0);
327 set_binary_interaction_double(i, j, "gammaT", 1.0);
328 set_binary_interaction_double(i, j, "betaV", 1.0);
329 set_binary_interaction_double(i, j, "gammaV", 1.0);
330 } else {
331 throw ValueError(format("mixing rule [%s] is not understood", model.c_str()));
332 }
333}
335void HelmholtzEOSMixtureBackend::set_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter,
336 const double value) {
337 // bound-check indices
338 if (i < 0 || i >= N) {
339 if (j < 0 || j >= N) {
340 throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N-1));
341 } else {
342 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N-1));
343 }
344 } else if (j < 0 || j >= N) {
345 throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, N-1));
346 }
347 if (parameter == "Fij") {
348 residual_helmholtz->Excess.F[i][j] = value;
349 residual_helmholtz->Excess.F[j][i] = value;
350 } else {
351 Reducing->set_binary_interaction_double(i, j, parameter, value);
352 }
354 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
355 it->get()->set_binary_interaction_double(i, j, parameter, value);
356 }
357};
359double HelmholtzEOSMixtureBackend::get_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter) {
360 // bound-check indices
361 if (i < 0 || i >= N) {
362 if (j < 0 || j >= N) {
363 throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N-1));
364 } else {
365 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N-1));
366 }
367 } else if (j < 0 || j >= N) {
368 throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, N-1));
369 }
370 if (parameter == "Fij") {
371 return residual_helmholtz->Excess.F[i][j];
372 } else {
373 return Reducing->get_binary_interaction_double(i, j, parameter);
374 }
375};
377//std::string HelmholtzEOSMixtureBackend::get_binary_interaction_string(const std::string &CAS1, const std::string &CAS2, const std::string &parameter){
378// return CoolProp::get_mixture_binary_pair_data(CAS1, CAS2, parameter);
379//}
381void HelmholtzEOSMixtureBackend::set_binary_interaction_string(const std::size_t i, const std::size_t j, const std::string& parameter,
382 const std::string& value) {
383 // bound-check indices
384 if (i < 0 || i >= N) {
385 if (j < 0 || j >= N) {
386 throw ValueError(format("Both indices i [%d] and j [%d] are out of bounds. Must be between 0 and %d.", i, j, N-1));
387 } else {
388 throw ValueError(format("Index i [%d] is out of bounds. Must be between 0 and %d.", i, N-1));
389 }
390 } else if (j < 0 || j >= N) {
391 throw ValueError(format("Index j [%d] is out of bounds. Must be between 0 and %d.", j, N-1));
392 }
393 if (parameter == "function") {
394 residual_helmholtz->Excess.DepartureFunctionMatrix[i][j].reset(get_departure_function(value));
395 residual_helmholtz->Excess.DepartureFunctionMatrix[j][i].reset(get_departure_function(value));
396 } else {
397 throw ValueError(format("Cannot process this string parameter [%s] in set_binary_interaction_string", parameter.c_str()));
398 }
400 for (std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>::iterator it = linked_states.begin(); it != linked_states.end(); ++it) {
401 it->get()->set_binary_interaction_string(i, j, parameter, value);
402 }
403};
404
405void HelmholtzEOSMixtureBackend::calc_change_EOS(const std::size_t i, const std::string& EOS_name) {
406
407 if (i < components.size()) {
408 CoolPropFluid& fluid = components[i];
409 EquationOfState& EOS = fluid.EOSVector[0];
410
411 if (EOS_name == "SRK" || EOS_name == "Peng-Robinson") {
412
413 // Get the parameters for the cubic EOS
414 CoolPropDbl Tc = EOS.reduce.T;
415 CoolPropDbl pc = EOS.reduce.p;
416 CoolPropDbl rhomolarc = EOS.reduce.rhomolar;
417 CoolPropDbl acentric = EOS.acentric;
418 CoolPropDbl R = 8.3144598;
419
420 // Remove the residual part
421 EOS.alphar.empty_the_EOS();
422 // Set the contribution
423 shared_ptr<AbstractCubic> ac;
424 if (EOS_name == "SRK") {
425 ac.reset(new SRK(Tc, pc, acentric, R));
426 } else {
427 ac.reset(new PengRobinson(Tc, pc, acentric, R));
428 }
429 ac->set_Tr(Tc);
430 ac->set_rhor(rhomolarc);
431 EOS.alphar.cubic = ResidualHelmholtzGeneralizedCubic(ac);
432 } else if (EOS_name == "XiangDeiters") {
433
434 // Get the parameters for the EOS
435 CoolPropDbl Tc = EOS.reduce.T;
436 CoolPropDbl pc = EOS.reduce.p;
437 CoolPropDbl rhomolarc = EOS.reduce.rhomolar;
438 CoolPropDbl acentric = EOS.acentric;
439 CoolPropDbl R = 8.3144598;
440
441 // Remove the residual part
442 EOS.alphar.empty_the_EOS();
443 // Set the Xiang & Deiters contribution
444 EOS.alphar.XiangDeiters = ResidualHelmholtzXiangDeiters(Tc, pc, rhomolarc, acentric, R);
445 }
446 } else {
447 throw ValueError(format("Index [%d] is invalid", i));
448 }
449 // Now do the same thing to the saturated liquid and vapor instances if possible
450 if (this->SatL) SatL->change_EOS(i, EOS_name);
451 if (this->SatV) SatV->change_EOS(i, EOS_name);
452}
453void HelmholtzEOSMixtureBackend::calc_phase_envelope(const std::string& type) {
454 // Clear the phase envelope data
455 PhaseEnvelope = PhaseEnvelopeData();
456 // Build the phase envelope
457 PhaseEnvelopeRoutines::build(*this, type);
458 // Finalize the phase envelope
459 PhaseEnvelopeRoutines::finalize(*this);
460};
461void HelmholtzEOSMixtureBackend::set_mixture_parameters() {
462 // Build the matrix of binary-pair reducing functions
463 MixtureParameters::set_mixture_parameters(*this);
464}
465void HelmholtzEOSMixtureBackend::update_states(void) {
466 CoolPropFluid& component = components[0];
467 EquationOfState& EOS = component.EOSVector[0];
468
469 // Clear the state class
470 clear();
471
472 // Calculate the new enthalpy and entropy values
473 update(DmolarT_INPUTS, EOS.hs_anchor.rhomolar, EOS.hs_anchor.T);
474 EOS.hs_anchor.hmolar = hmolar();
475 EOS.hs_anchor.smolar = smolar();
476
477 // Calculate the new enthalpy and entropy values at the reducing state
478 update(DmolarT_INPUTS, EOS.reduce.rhomolar, EOS.reduce.T);
479 EOS.reduce.hmolar = hmolar();
480 EOS.reduce.smolar = smolar();
481
482 // Clear again just to be sure
483 clear();
484}
485const CoolProp::SimpleState& HelmholtzEOSMixtureBackend::calc_state(const std::string& state) {
486 if (is_pure_or_pseudopure) {
487 if (!state.compare("hs_anchor")) {
488 return components[0].EOS().hs_anchor;
489 } else if (!state.compare("max_sat_T")) {
490 return components[0].EOS().max_sat_T;
491 } else if (!state.compare("max_sat_p")) {
492 return components[0].EOS().max_sat_p;
493 } else if (!state.compare("reducing")) {
494 return components[0].EOS().reduce;
495 } else if (!state.compare("critical")) {
496 return components[0].crit;
497 } else if (!state.compare("triple_liquid")) {
498 return components[0].triple_liquid;
499 } else if (!state.compare("triple_vapor")) {
500 return components[0].triple_vapor;
501 } else {
502 throw ValueError(format("This state [%s] is invalid to calc_state", state.c_str()));
503 }
504 } else {
505 if (!state.compare("critical")) {
506 return _critical;
507 } else {
508 throw ValueError(format("calc_state not supported for mixtures"));
509 }
510 }
511};
512CoolPropDbl HelmholtzEOSMixtureBackend::calc_acentric_factor(void) {
513 if (is_pure_or_pseudopure) {
514 return components[0].EOS().acentric;
515 } else {
516 throw ValueError("acentric factor cannot be calculated for mixtures");
517 }
518}
519CoolPropDbl HelmholtzEOSMixtureBackend::calc_gas_constant(void) {
520 if (is_pure_or_pseudopure) {
521 return components[0].gas_constant();
522 } else {
523 if (get_config_bool(NORMALIZE_GAS_CONSTANTS)) {
524 return get_config_double(R_U_CODATA);
525 } else {
526 // mass fraction weighted average of the components
527 double summer = 0;
528 for (unsigned int i = 0; i < components.size(); ++i) {
529 summer += mole_fractions[i] * components[i].gas_constant();
530 }
531 return summer;
532 }
533 }
534}
535CoolPropDbl HelmholtzEOSMixtureBackend::calc_molar_mass(void) {
536 double summer = 0;
537 for (unsigned int i = 0; i < components.size(); ++i) {
538 summer += mole_fractions[i] * components[i].molar_mass();
539 }
540 return summer;
541}
542CoolPropDbl HelmholtzEOSMixtureBackend::calc_saturation_ancillary(parameters param, int Q, parameters given, double value) {
543 if (is_pure_or_pseudopure) {
544 if (param == iP && given == iT) {
545 // p = f(T), direct evaluation
546 switch (Q) {
547 case 0:
548 return components[0].ancillaries.pL.evaluate(value);
549 case 1:
550 return components[0].ancillaries.pV.evaluate(value);
551 }
552 } else if (param == iT && given == iP) {
553 // T = f(p), inverse evaluation
554 switch (Q) {
555 case 0:
556 return components[0].ancillaries.pL.invert(value);
557 case 1:
558 return components[0].ancillaries.pV.invert(value);
559 }
560 } else if (param == iDmolar && given == iT) {
561 // rho = f(T), inverse evaluation
562 switch (Q) {
563 case 0:
564 return components[0].ancillaries.rhoL.evaluate(value);
565 case 1:
566 return components[0].ancillaries.rhoV.evaluate(value);
567 }
568 } else if (param == iT && given == iDmolar) {
569 // T = f(rho), inverse evaluation
570 switch (Q) {
571 case 0:
572 return components[0].ancillaries.rhoL.invert(value);
573 case 1:
574 return components[0].ancillaries.rhoV.invert(value);
575 }
576 } else if (param == isurface_tension && given == iT) {
577 return components[0].ancillaries.surface_tension.evaluate(value);
578 } else {
579 throw ValueError(format("calc of %s given %s is invalid in calc_saturation_ancillary", get_parameter_information(param, "short").c_str(),
580 get_parameter_information(given, "short").c_str()));
581 }
582
583 throw ValueError(format("Q [%d] is invalid in calc_saturation_ancillary", Q));
584 } else {
585 throw NotImplementedError(format("calc_saturation_ancillary not implemented for mixtures"));
586 }
587}
588
589CoolPropDbl HelmholtzEOSMixtureBackend::calc_melting_line(int param, int given, CoolPropDbl value) {
590 if (is_pure_or_pseudopure) {
591 return components[0].ancillaries.melting_line.evaluate(param, given, value);
592 } else {
593 throw NotImplementedError(format("calc_melting_line not implemented for mixtures"));
594 }
595}
596CoolPropDbl HelmholtzEOSMixtureBackend::calc_surface_tension(void) {
597 if (is_pure_or_pseudopure) {
598 if ((_phase == iphase_twophase) || (_phase == iphase_critical_point)) { // if within the two phase region or at critical point
599 return components[0].ancillaries.surface_tension.evaluate(T()); // calculate surface tension and return
600 } else { // else state point not in the two phase region
601 throw ValueError(format("surface tension is only defined within the two-phase region; Try PQ or QT inputs")); // throw error
602 }
603 } else {
604 throw NotImplementedError(format("surface tension not implemented for mixtures"));
605 }
606}
607CoolPropDbl HelmholtzEOSMixtureBackend::calc_viscosity_dilute(void) {
608 if (is_pure_or_pseudopure) {
609 CoolPropDbl eta_dilute;
610 switch (components[0].transport.viscosity_dilute.type) {
611 case ViscosityDiluteVariables::VISCOSITY_DILUTE_KINETIC_THEORY:
612 eta_dilute = TransportRoutines::viscosity_dilute_kinetic_theory(*this);
613 break;
614 case ViscosityDiluteVariables::VISCOSITY_DILUTE_COLLISION_INTEGRAL:
615 eta_dilute = TransportRoutines::viscosity_dilute_collision_integral(*this);
616 break;
617 case ViscosityDiluteVariables::VISCOSITY_DILUTE_POWERS_OF_T:
618 eta_dilute = TransportRoutines::viscosity_dilute_powers_of_T(*this);
619 break;
620 case ViscosityDiluteVariables::VISCOSITY_DILUTE_POWERS_OF_TR:
621 eta_dilute = TransportRoutines::viscosity_dilute_powers_of_Tr(*this);
622 break;
623 case ViscosityDiluteVariables::VISCOSITY_DILUTE_COLLISION_INTEGRAL_POWERS_OF_TSTAR:
624 eta_dilute = TransportRoutines::viscosity_dilute_collision_integral_powers_of_T(*this);
625 break;
626 case ViscosityDiluteVariables::VISCOSITY_DILUTE_ETHANE:
627 eta_dilute = TransportRoutines::viscosity_dilute_ethane(*this);
628 break;
629 case ViscosityDiluteVariables::VISCOSITY_DILUTE_CYCLOHEXANE:
630 eta_dilute = TransportRoutines::viscosity_dilute_cyclohexane(*this);
631 break;
632 case ViscosityDiluteVariables::VISCOSITY_DILUTE_CO2_LAESECKE_JPCRD_2017:
633 eta_dilute = TransportRoutines::viscosity_dilute_CO2_LaeseckeJPCRD2017(*this);
634 break;
635 default:
636 throw ValueError(
637 format("dilute viscosity type [%d] is invalid for fluid %s", components[0].transport.viscosity_dilute.type, name().c_str()));
638 }
639 return eta_dilute;
640 } else {
641 throw NotImplementedError(format("dilute viscosity not implemented for mixtures"));
642 }
643}
644CoolPropDbl HelmholtzEOSMixtureBackend::calc_viscosity_background() {
645 CoolPropDbl eta_dilute = calc_viscosity_dilute(), initial_density = 0, residual = 0;
646 return calc_viscosity_background(eta_dilute, initial_density, residual);
647}
648CoolPropDbl HelmholtzEOSMixtureBackend::calc_viscosity_background(CoolPropDbl eta_dilute, CoolPropDbl& initial_density, CoolPropDbl& residual) {
649
650 switch (components[0].transport.viscosity_initial.type) {
651 case ViscosityInitialDensityVariables::VISCOSITY_INITIAL_DENSITY_RAINWATER_FRIEND: {
652 CoolPropDbl B_eta_initial = TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend(*this);
653 CoolPropDbl rho = rhomolar();
654 initial_density = eta_dilute * B_eta_initial * rho; //TODO: Check units once AMTG
655 break;
656 }
657 case ViscosityInitialDensityVariables::VISCOSITY_INITIAL_DENSITY_EMPIRICAL: {
658 initial_density = TransportRoutines::viscosity_initial_density_dependence_empirical(*this);
659 break;
660 }
661 case ViscosityInitialDensityVariables::VISCOSITY_INITIAL_DENSITY_NOT_SET: {
662 break;
663 }
664 }
665
666 // Higher order terms
667 switch (components[0].transport.viscosity_higher_order.type) {
668 case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_BATSCHINKI_HILDEBRAND:
669 residual = TransportRoutines::viscosity_higher_order_modified_Batschinski_Hildebrand(*this);
670 break;
671 case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_FRICTION_THEORY:
672 residual = TransportRoutines::viscosity_higher_order_friction_theory(*this);
673 break;
674 case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_HYDROGEN:
675 residual = TransportRoutines::viscosity_hydrogen_higher_order_hardcoded(*this);
676 break;
677 case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_TOLUENE:
678 residual = TransportRoutines::viscosity_toluene_higher_order_hardcoded(*this);
679 break;
680 case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_HEXANE:
681 residual = TransportRoutines::viscosity_hexane_higher_order_hardcoded(*this);
682 break;
683 case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_HEPTANE:
684 residual = TransportRoutines::viscosity_heptane_higher_order_hardcoded(*this);
685 break;
686 case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_ETHANE:
687 residual = TransportRoutines::viscosity_ethane_higher_order_hardcoded(*this);
688 break;
689 case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_BENZENE:
690 residual = TransportRoutines::viscosity_benzene_higher_order_hardcoded(*this);
691 break;
692 case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_CO2_LAESECKE_JPCRD_2017:
693 residual = TransportRoutines::viscosity_CO2_higher_order_hardcoded_LaeseckeJPCRD2017(*this);
694 break;
695 default:
696 throw ValueError(
697 format("higher order viscosity type [%d] is invalid for fluid %s", components[0].transport.viscosity_dilute.type, name().c_str()));
698 }
699
700 return initial_density + residual;
701}
702
703CoolPropDbl HelmholtzEOSMixtureBackend::calc_viscosity(void) {
704 if (is_pure_or_pseudopure) {
705 CoolPropDbl dilute = 0, initial_density = 0, residual = 0, critical = 0;
706 calc_viscosity_contributions(dilute, initial_density, residual, critical);
707 return dilute + initial_density + residual + critical;
708 } else {
709 set_warning_string("Mixture model for viscosity is highly approximate");
710 CoolPropDbl summer = 0;
711 for (std::size_t i = 0; i < mole_fractions.size(); ++i) {
712 shared_ptr<HelmholtzEOSBackend> HEOS(new HelmholtzEOSBackend(components[i]));
713 HEOS->update(DmolarT_INPUTS, _rhomolar, _T);
714 summer += mole_fractions[i] * log(HEOS->viscosity());
715 }
716 return exp(summer);
717 }
718}
719void HelmholtzEOSMixtureBackend::calc_viscosity_contributions(CoolPropDbl& dilute, CoolPropDbl& initial_density, CoolPropDbl& residual,
720 CoolPropDbl& critical) {
721 if (is_pure_or_pseudopure) {
722 // Reset the variables
723 dilute = 0;
724 initial_density = 0;
725 residual = 0;
726 critical = 0;
727
728 // Get a reference for code cleanness
729 CoolPropFluid& component = components[0];
730
731 if (!component.transport.viscosity_model_provided) {
732 throw ValueError(format("Viscosity model is not available for this fluid"));
733 }
734
735 // Check if using ECS
736 if (component.transport.viscosity_using_ECS) {
737 // Get reference fluid name
738 std::string fluid_name = component.transport.viscosity_ecs.reference_fluid;
739 std::vector<std::string> names(1, fluid_name);
740 // Get a managed pointer to the reference fluid for ECS
741 shared_ptr<HelmholtzEOSMixtureBackend> ref_fluid(new HelmholtzEOSMixtureBackend(names));
742 // Get the viscosity using ECS and stick in the critical value
743 critical = TransportRoutines::viscosity_ECS(*this, *ref_fluid);
744 return;
745 }
746
747 // Check if using Chung model
748 if (component.transport.viscosity_using_Chung) {
749 // Get the viscosity using ECS and stick in the critical value
750 critical = TransportRoutines::viscosity_Chung(*this);
751 return;
752 }
753
754 // Check if using rho*sr model
755 if (component.transport.viscosity_using_rhosr) {
756 // Get the viscosity using rho*sr model and stick in the critical value
757 critical = TransportRoutines::viscosity_rhosr(*this);
758 return;
759 }
760
761 if (component.transport.hardcoded_viscosity != CoolProp::TransportPropertyData::VISCOSITY_NOT_HARDCODED) {
762 // Evaluate hardcoded model and stick in the critical value
763 switch (component.transport.hardcoded_viscosity) {
764 case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_WATER:
765 critical = TransportRoutines::viscosity_water_hardcoded(*this);
766 break;
767 case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_HEAVYWATER:
768 critical = TransportRoutines::viscosity_heavywater_hardcoded(*this);
769 break;
770 case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_HELIUM:
771 critical = TransportRoutines::viscosity_helium_hardcoded(*this);
772 break;
773 case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_R23:
774 critical = TransportRoutines::viscosity_R23_hardcoded(*this);
775 break;
776 case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_METHANOL:
777 critical = TransportRoutines::viscosity_methanol_hardcoded(*this);
778 break;
779 case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_M_XYLENE:
780 critical = TransportRoutines::viscosity_m_xylene_hardcoded(*this);
781 break;
782 case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_O_XYLENE:
783 critical = TransportRoutines::viscosity_o_xylene_hardcoded(*this);
784 break;
785 case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_P_XYLENE:
786 critical = TransportRoutines::viscosity_p_xylene_hardcoded(*this);
787 break;
788 default:
789 throw ValueError(
790 format("hardcoded viscosity type [%d] is invalid for fluid %s", component.transport.hardcoded_viscosity, name().c_str()));
791 }
792 return;
793 }
794
795 // -------------------------
796 // Normal evaluation
797 // -------------------------
798
799 // Dilute part
800 dilute = calc_viscosity_dilute();
801
802 // Background viscosity given by the sum of the initial density dependence and higher order terms
803 calc_viscosity_background(dilute, initial_density, residual);
804
805 // Critical part (no fluids have critical enhancement for viscosity currently)
806 critical = 0;
807 } else {
808 throw ValueError("calc_viscosity_contributions invalid for mixtures");
809 }
810}
811void HelmholtzEOSMixtureBackend::calc_conductivity_contributions(CoolPropDbl& dilute, CoolPropDbl& initial_density, CoolPropDbl& residual,
812 CoolPropDbl& critical) {
813 if (is_pure_or_pseudopure) {
814 // Reset the variables
815 dilute = 0;
816 initial_density = 0;
817 residual = 0;
818 critical = 0;
819
820 // Get a reference for code cleanness
821 CoolPropFluid& component = components[0];
822
823 if (!component.transport.conductivity_model_provided) {
824 throw ValueError(format("Thermal conductivity model is not available for this fluid"));
825 }
826
827 // Check if using ECS
828 if (component.transport.conductivity_using_ECS) {
829 // Get reference fluid name
830 std::string fluid_name = component.transport.conductivity_ecs.reference_fluid;
831 std::vector<std::string> name(1, fluid_name);
832 // Get a managed pointer to the reference fluid for ECS
833 shared_ptr<HelmholtzEOSMixtureBackend> ref_fluid(new HelmholtzEOSMixtureBackend(name));
834 // Get the viscosity using ECS and store in initial_density (not normally used);
835 initial_density = TransportRoutines::conductivity_ECS(*this, *ref_fluid); // Warning: not actually initial_density
836 return;
837 }
838
839 if (component.transport.hardcoded_conductivity != CoolProp::TransportPropertyData::CONDUCTIVITY_NOT_HARDCODED) {
840 // Evaluate hardcoded model and deposit in initial_density variable
841 // Warning: not actually initial_density
842 switch (component.transport.hardcoded_conductivity) {
843 case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_WATER:
844 initial_density = TransportRoutines::conductivity_hardcoded_water(*this);
845 break;
846 case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_HEAVYWATER:
847 initial_density = TransportRoutines::conductivity_hardcoded_heavywater(*this);
848 break;
849 case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_R23:
850 initial_density = TransportRoutines::conductivity_hardcoded_R23(*this);
851 break;
852 case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_HELIUM:
853 initial_density = TransportRoutines::conductivity_hardcoded_helium(*this);
854 break;
855 case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_METHANE:
856 initial_density = TransportRoutines::conductivity_hardcoded_methane(*this);
857 break;
858 default:
859 throw ValueError(format("hardcoded conductivity type [%d] is invalid for fluid %s",
860 components[0].transport.hardcoded_conductivity, name().c_str()));
861 }
862 return;
863 }
864
865 // -------------------------
866 // Normal evaluation
867 // -------------------------
868
869 // Dilute part
870 switch (component.transport.conductivity_dilute.type) {
871 case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_RATIO_POLYNOMIALS:
872 dilute = TransportRoutines::conductivity_dilute_ratio_polynomials(*this);
873 break;
874 case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_ETA0_AND_POLY:
875 dilute = TransportRoutines::conductivity_dilute_eta0_and_poly(*this);
876 break;
877 case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_CO2:
878 dilute = TransportRoutines::conductivity_dilute_hardcoded_CO2(*this);
879 break;
880 case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_CO2_HUBER_JPCRD_2016:
881 dilute = TransportRoutines::conductivity_dilute_hardcoded_CO2_HuberJPCRD2016(*this);
882 break;
883 case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_ETHANE:
884 dilute = TransportRoutines::conductivity_dilute_hardcoded_ethane(*this);
885 break;
886 case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_NONE:
887 dilute = 0.0;
888 break;
889 default:
890 throw ValueError(
891 format("dilute conductivity type [%d] is invalid for fluid %s", components[0].transport.conductivity_dilute.type, name().c_str()));
892 }
893
894 // Residual part
895 residual = calc_conductivity_background();
896
897 // Critical part
898 switch (component.transport.conductivity_critical.type) {
899 case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_SIMPLIFIED_OLCHOWY_SENGERS:
900 critical = TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(*this);
901 break;
902 case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_R123:
903 critical = TransportRoutines::conductivity_critical_hardcoded_R123(*this);
904 break;
905 case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_AMMONIA:
906 critical = TransportRoutines::conductivity_critical_hardcoded_ammonia(*this);
907 break;
908 case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_NONE:
909 critical = 0.0;
910 break;
911 case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_CARBONDIOXIDE_SCALABRIN_JPCRD_2006:
912 critical = TransportRoutines::conductivity_critical_hardcoded_CO2_ScalabrinJPCRD2006(*this);
913 break;
914 default:
915 throw ValueError(
916 format("critical conductivity type [%d] is invalid for fluid %s", components[0].transport.viscosity_dilute.type, name().c_str()));
917 }
918 } else {
919 throw ValueError("calc_conductivity_contributions invalid for mixtures");
920 }
921};
922
923CoolPropDbl HelmholtzEOSMixtureBackend::calc_conductivity_background(void) {
924 // Residual part
925 CoolPropDbl lambda_residual = _HUGE;
926 switch (components[0].transport.conductivity_residual.type) {
927 case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL:
928 lambda_residual = TransportRoutines::conductivity_residual_polynomial(*this);
929 break;
930 case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL_AND_EXPONENTIAL:
931 lambda_residual = TransportRoutines::conductivity_residual_polynomial_and_exponential(*this);
932 break;
933 default:
934 throw ValueError(
935 format("residual conductivity type [%d] is invalid for fluid %s", components[0].transport.conductivity_residual.type, name().c_str()));
936 }
937 return lambda_residual;
938}
939CoolPropDbl HelmholtzEOSMixtureBackend::calc_conductivity(void) {
940 if (is_pure_or_pseudopure) {
941 CoolPropDbl dilute = 0, initial_density = 0, residual = 0, critical = 0;
942 calc_conductivity_contributions(dilute, initial_density, residual, critical);
943 return dilute + initial_density + residual + critical;
944 } else {
945 set_warning_string("Mixture model for conductivity is highly approximate");
946 CoolPropDbl summer = 0;
947 for (std::size_t i = 0; i < mole_fractions.size(); ++i) {
948 shared_ptr<HelmholtzEOSBackend> HEOS(new HelmholtzEOSBackend(components[i]));
949 HEOS->update(DmolarT_INPUTS, _rhomolar, _T);
950 summer += mole_fractions[i] * HEOS->conductivity();
951 }
952 return summer;
953 }
954}
955void HelmholtzEOSMixtureBackend::calc_conformal_state(const std::string& reference_fluid, CoolPropDbl& T, CoolPropDbl& rhomolar) {
956 shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> REF(new CoolProp::HelmholtzEOSBackend(reference_fluid));
957
958 if (T < 0 && rhomolar < 0) {
959 // Collect some parameters
960 CoolPropDbl Tc = T_critical(), Tc0 = REF->T_critical(), rhocmolar = rhomolar_critical(), rhocmolar0 = REF->rhomolar_critical();
961
962 // Starting guess values for shape factors
963 CoolPropDbl theta = 1;
964 CoolPropDbl phi = 1;
965
966 // The equivalent substance reducing ratios
967 CoolPropDbl f = Tc / Tc0 * theta;
968 CoolPropDbl h = rhocmolar0 / rhocmolar * phi; // Must be the ratio of MOLAR densities!!
969
970 // Starting guesses for conformal state
971 T = this->T() / f;
972 rhomolar = this->rhomolar() * h;
973 }
974
975 TransportRoutines::conformal_state_solver(*this, *REF, T, rhomolar);
976}
977CoolPropDbl HelmholtzEOSMixtureBackend::calc_Ttriple(void) {
978 double summer = 0;
979 for (unsigned int i = 0; i < components.size(); ++i) {
980 summer += mole_fractions[i] * components[i].EOS().Ttriple;
981 }
982 return summer;
983}
984CoolPropDbl HelmholtzEOSMixtureBackend::calc_p_triple(void) {
985 double summer = 0;
986 for (unsigned int i = 0; i < components.size(); ++i) {
987 summer += mole_fractions[i] * components[i].EOS().ptriple;
988 }
989 return summer;
990}
991std::string HelmholtzEOSMixtureBackend::calc_name(void) {
992 if (components.size() != 1) {
993 throw ValueError(format("calc_name is only valid for pure and pseudo-pure fluids, %d components", components.size()));
994 } else {
995 return components[0].name;
996 }
997}
998
999CoolPropDbl HelmholtzEOSMixtureBackend::calc_saturated_liquid_keyed_output(parameters key) {
1000 if ((key == iDmolar) && _rhoLmolar) return _rhoLmolar;
1001 if (!SatL) throw ValueError("The saturated liquid state has not been set.");
1002 return SatL->keyed_output(key);
1003}
1004CoolPropDbl HelmholtzEOSMixtureBackend::calc_saturated_vapor_keyed_output(parameters key) {
1005 if ((key == iDmolar) && _rhoVmolar) return _rhoVmolar;
1006 if (!SatV) throw ValueError("The saturated vapor state has not been set.");
1007 return SatV->keyed_output(key);
1008}
1009
1010void HelmholtzEOSMixtureBackend::calc_ideal_curve(const std::string& type, std::vector<double>& T, std::vector<double>& p) {
1011 if (type == "Joule-Thomson") {
1012 JouleThomsonCurveTracer JTCT(this, 1e5, 800);
1013 JTCT.trace(T, p);
1014 } else if (type == "Joule-Inversion") {
1015 JouleInversionCurveTracer JICT(this, 1e5, 800);
1016 JICT.trace(T, p);
1017 } else if (type == "Ideal") {
1018 IdealCurveTracer ICT(this, 1e5, 800);
1019 ICT.trace(T, p);
1020 } else if (type == "Boyle") {
1021 BoyleCurveTracer BCT(this, 1e5, 800);
1022 BCT.trace(T, p);
1023 } else {
1024 throw ValueError(format("Invalid ideal curve type: %s", type.c_str()));
1025 }
1026};
1027std::vector<std::string> HelmholtzEOSMixtureBackend::calc_fluid_names(void) {
1028 std::vector<std::string> out;
1029 for (std::size_t i = 0; i < components.size(); ++i) {
1030 out.push_back(components[i].name);
1031 }
1032 return out;
1033}
1034CoolPropDbl HelmholtzEOSMixtureBackend::calc_ODP(void) {
1035 if (components.size() != 1) {
1036 throw ValueError(format("For now, calc_ODP is only valid for pure and pseudo-pure fluids, %d components", components.size()));
1037 } else {
1038 CoolPropDbl v = components[0].environment.ODP;
1039 if (!ValidNumber(v) || v < 0) {
1040 throw ValueError(format("ODP value is not specified or invalid"));
1041 }
1042 return v;
1043 }
1044}
1045CoolPropDbl HelmholtzEOSMixtureBackend::calc_GWP20(void) {
1046 if (components.size() != 1) {
1047 throw ValueError(format("For now, calc_GWP20 is only valid for pure and pseudo-pure fluids, %d components", components.size()));
1048 } else {
1049 CoolPropDbl v = components[0].environment.GWP20;
1050 if (!ValidNumber(v) || v < 0) {
1051 throw ValueError(format("GWP20 value is not specified or invalid"));
1052 }
1053 return v;
1054 }
1055}
1056CoolPropDbl HelmholtzEOSMixtureBackend::calc_GWP100(void) {
1057 if (components.size() != 1) {
1058 throw ValueError(format("For now, calc_GWP100 is only valid for pure and pseudo-pure fluids, %d components", components.size()));
1059 } else {
1060 CoolPropDbl v = components[0].environment.GWP100;
1061 if (!ValidNumber(v) || v < 0) {
1062 throw ValueError(format("GWP100 value is not specified or invalid"));
1063 }
1064 return v;
1065 }
1066}
1067CoolPropDbl HelmholtzEOSMixtureBackend::calc_GWP500(void) {
1068 if (components.size() != 1) {
1069 throw ValueError(format("For now, calc_GWP500 is only valid for pure and pseudo-pure fluids, %d components", components.size()));
1070 } else {
1071 CoolPropDbl v = components[0].environment.GWP500;
1072 if (!ValidNumber(v) || v < 0) {
1073 throw ValueError(format("GWP500 value is not specified or invalid"));
1074 }
1075 return v;
1076 }
1077}
1078CoolPropDbl HelmholtzEOSMixtureBackend::calc_T_critical(void) {
1079 if (components.size() != 1) {
1080 std::vector<CriticalState> critpts = calc_all_critical_points();
1081 if (critpts.size() == 1) {
1082 //if (!critpts[0].stable){ throw ValueError(format("found one critical point but critical point is not stable")); }
1083 return critpts[0].T;
1084 } else {
1085 throw ValueError(format("critical point finding routine found %d critical points", critpts.size()));
1086 }
1087 } else {
1088 if (get_config_bool(ENABLE_SUPERANCILLARIES) && is_pure()){
1089 auto& optsuperanc = get_superanc_optional();
1090 if (optsuperanc){
1091 return optsuperanc.value().get_Tcrit_num(); // from the numerical critical point satisfying dp/drho|T = d2p/drho2|T = 0
1092 }
1093 }
1094 return components[0].crit.T;
1095 }
1096}
1097CoolPropDbl HelmholtzEOSMixtureBackend::calc_p_critical(void) {
1098 if (components.size() != 1) {
1099 std::vector<CriticalState> critpts = calc_all_critical_points();
1100 if (critpts.size() == 1) {
1101 //if (!critpts[0].stable){ throw ValueError(format("found one critical point but critical point is not stable")); }
1102 return critpts[0].p;
1103 } else {
1104 throw ValueError(format("critical point finding routine found %d critical points", critpts.size()));
1105 }
1106 } else {
1107 if (get_config_bool(ENABLE_SUPERANCILLARIES) && is_pure()){
1108 auto& optsuperanc = get_superanc_optional();
1109 if (optsuperanc){
1110 return optsuperanc.value().get_pmax(); // from the numerical critical point satisfying dp/drho|T = d2p/drho2|T = 0
1111 }
1112 }
1113 return components[0].crit.p;
1114 }
1115}
1116CoolPropDbl HelmholtzEOSMixtureBackend::calc_rhomolar_critical(void) {
1117 if (components.size() != 1) {
1118 std::vector<CriticalState> critpts = calc_all_critical_points();
1119 if (critpts.size() == 1) {
1120 //if (!critpts[0].stable){ throw ValueError(format("found one critical point but critical point is not stable")); }
1121 return critpts[0].rhomolar;
1122 } else {
1123 throw ValueError(format("critical point finding routine found %d critical points", critpts.size()));
1124 }
1125 } else {
1126 if (get_config_bool(ENABLE_SUPERANCILLARIES) && is_pure()){
1127 auto& optsuperanc = get_superanc_optional();
1128 if (optsuperanc){
1129 return optsuperanc.value().get_rhocrit_num(); // from the numerical critical point satisfying dp/drho|T = d2p/drho2|T = 0
1130 }
1131 }
1132 return components[0].crit.rhomolar;
1133 }
1134}
1135CoolPropDbl HelmholtzEOSMixtureBackend::calc_pmax_sat(void) {
1136 if (is_pure_or_pseudopure) {
1137 if (components[0].EOS().pseudo_pure) {
1138 return components[0].EOS().max_sat_p.p;
1139 } else {
1140 return p_critical();
1141 }
1142 } else {
1143 throw ValueError("calc_pmax_sat not yet defined for mixtures");
1144 }
1145}
1146CoolPropDbl HelmholtzEOSMixtureBackend::calc_Tmax_sat(void) {
1147 if (is_pure_or_pseudopure) {
1148 if (components[0].EOS().pseudo_pure) {
1149 double Tmax_sat = components[0].EOS().max_sat_T.T;
1150 if (!ValidNumber(Tmax_sat)) {
1151 return T_critical();
1152 } else {
1153 return Tmax_sat;
1154 }
1155 } else {
1156 return T_critical();
1157 }
1158 } else {
1159 throw ValueError("calc_Tmax_sat not yet defined for mixtures");
1160 }
1161}
1162
1163void HelmholtzEOSMixtureBackend::calc_Tmin_sat(CoolPropDbl& Tmin_satL, CoolPropDbl& Tmin_satV) {
1164 if (is_pure_or_pseudopure) {
1165 Tmin_satL = components[0].EOS().sat_min_liquid.T;
1166 Tmin_satV = components[0].EOS().sat_min_vapor.T;
1167 return;
1168 } else {
1169 throw ValueError("calc_Tmin_sat not yet defined for mixtures");
1170 }
1171}
1172
1173void HelmholtzEOSMixtureBackend::calc_pmin_sat(CoolPropDbl& pmin_satL, CoolPropDbl& pmin_satV) {
1174 if (is_pure_or_pseudopure) {
1175 pmin_satL = components[0].EOS().sat_min_liquid.p;
1176 pmin_satV = components[0].EOS().sat_min_vapor.p;
1177 return;
1178 } else {
1179 throw ValueError("calc_pmin_sat not yet defined for mixtures");
1180 }
1181}
1182
1183// Minimum allowed saturation temperature the maximum of the saturation temperatures of liquid and vapor
1184// For pure fluids, both values are the same, for pseudo-pure they are probably the same, for mixtures they are definitely not the same
1185
1186CoolPropDbl HelmholtzEOSMixtureBackend::calc_Tmax(void) {
1187 double summer = 0;
1188 for (unsigned int i = 0; i < components.size(); ++i) {
1189 summer += mole_fractions[i] * components[i].EOS().limits.Tmax;
1190 }
1191 return summer;
1192}
1193CoolPropDbl HelmholtzEOSMixtureBackend::calc_Tmin(void) {
1194 double summer = 0;
1195 for (unsigned int i = 0; i < components.size(); ++i) {
1196 summer += mole_fractions[i] * components[i].EOS().limits.Tmin;
1197 }
1198 return summer;
1199}
1200CoolPropDbl HelmholtzEOSMixtureBackend::calc_pmax(void) {
1201 double summer = 0;
1202 for (unsigned int i = 0; i < components.size(); ++i) {
1203 summer += mole_fractions[i] * components[i].EOS().limits.pmax;
1204 }
1205 return summer;
1206}
1207
1208void HelmholtzEOSMixtureBackend::update_QT_pure_superanc(CoolPropDbl Q, CoolPropDbl T) {
1209 clear();
1210 _Q = Q;
1211 _T = T;
1212 if (!is_pure()){
1213 throw ValueError(format("Is not a pure fluid"));
1214 }
1215 if (!get_config_bool(ENABLE_SUPERANCILLARIES)){
1216 throw ValueError(format("Superancillaries are not enabled"));
1217 }
1218 auto& optsuperanc = get_superanc_optional();
1219 if (!optsuperanc){
1220 throw ValueError(format("Superancillaries not available for this fluid"));
1221 }
1222
1223 auto& superanc = optsuperanc.value();
1224 CoolPropDbl Tcrit_num = superanc.get_Tcrit_num();
1225 if (T > Tcrit_num){
1226 throw ValueError(format("Temperature to QT_flash [%0.8Lg K] may not be above the numerical critical point of %0.15Lg K", T, Tcrit_num));
1227 }
1228 auto rhoL = superanc.eval_sat(T, 'D', 0);
1229 auto rhoV = superanc.eval_sat(T, 'D', 1);
1230 auto p = superanc.eval_sat(T, 'P', 1);
1231 SatL->update_TDmolarP_unchecked(T, rhoL, p);
1232 SatV->update_TDmolarP_unchecked(T, rhoV, p);
1233 _p = p;
1234 _rhomolar = 1 / (Q / rhoV + (1 - Q) / rhoL);
1235 _phase = iphase_twophase;
1236
1237 post_update(false /*optional_checks*/);
1238}
1239
1240
1241void HelmholtzEOSMixtureBackend::update_TDmolarP_unchecked(CoolPropDbl T, CoolPropDbl rhomolar, CoolPropDbl p) {
1242
1243 clear();
1244
1245 _rhomolar = rhomolar;
1246 _T = T;
1247 _p = p;
1248
1249 // Cleanup
1250 bool optional_checks = false;
1251 post_update(optional_checks);
1252}
1253
1254
1255
1256void HelmholtzEOSMixtureBackend::update_DmolarT_direct(CoolPropDbl rhomolar, CoolPropDbl T) {
1257 // TODO: This is just a quick fix for #878 - should be done more systematically
1258 const CoolPropDbl rhomolar_min = 0;
1259 const CoolPropDbl T_min = 0;
1260
1261 if (rhomolar < rhomolar_min) {
1262 throw ValueError(format("The molar density of %f mol/m3 is below the minimum of %f mol/m3", rhomolar, rhomolar_min));
1263 }
1264
1265 if (T < T_min) {
1266 throw ValueError(format("The temperature of %f K is below the minimum of %f K", T, T_min));
1267 }
1268
1269 CoolProp::input_pairs pair = DmolarT_INPUTS;
1270 // Set up the state
1271 pre_update(pair, rhomolar, T);
1272
1273 _rhomolar = rhomolar;
1274 _T = T;
1275 _p = calc_pressure();
1276
1277 // Cleanup
1278 bool optional_checks = false;
1279 post_update(optional_checks);
1280}
1281
1282void HelmholtzEOSMixtureBackend::update_HmolarQ_with_guessT(CoolPropDbl hmolar, CoolPropDbl Q, CoolPropDbl Tguess) {
1283 CoolProp::input_pairs pair = CoolProp::HmolarQ_INPUTS;
1284 // Set up the state
1285 pre_update(pair, hmolar, Q);
1286
1287 _hmolar = hmolar;
1288 _Q = Q;
1289 FlashRoutines::HQ_flash(*this, Tguess);
1290
1291 // Cleanup
1292 post_update();
1293}
1294void HelmholtzEOSMixtureBackend::update_internal(HelmholtzEOSMixtureBackend& HEOS) {
1295 this->_hmolar = HEOS.hmolar();
1296 this->_smolar = HEOS.smolar();
1297 this->_T = HEOS.T();
1298 this->_umolar = HEOS.umolar();
1299 this->_p = HEOS.p();
1300 this->_rhomolar = HEOS.rhomolar();
1301 this->_Q = HEOS.Q();
1302 this->_phase = HEOS.phase();
1303
1304 // Copy the derivatives as well
1305}
1306void HelmholtzEOSMixtureBackend::update_TP_guessrho(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomolar_guess) {
1307 CoolProp::input_pairs pair = PT_INPUTS;
1308 // Set up the state
1309 pre_update(pair, p, T);
1310
1311 // Do the flash call to find rho = f(T,p)
1312 CoolPropDbl rhomolar = solver_rho_Tp(T, p, rhomolar_guess);
1313
1314 // Update the class with the new calculated density
1315 update_DmolarT_direct(rhomolar, T);
1316
1317 // Skip the cleanup, already done in update_DmolarT_direct
1318}
1319
1320void HelmholtzEOSMixtureBackend::pre_update(CoolProp::input_pairs& input_pair, CoolPropDbl& value1, CoolPropDbl& value2) {
1321 // Clear the state
1322 clear();
1323
1324 if (is_pure_or_pseudopure == false && mole_fractions.size() == 0) {
1325 throw ValueError("Mole fractions must be set");
1326 }
1327
1328 // If the inputs are in mass units, convert them to molar units
1329 mass_to_molar_inputs(input_pair, value1, value2);
1330
1331 // Set the mole-fraction weighted gas constant for the mixture
1332 // (or the pure/pseudo-pure fluid) if it hasn't been set yet
1333 gas_constant();
1334
1335 // Calculate and cache the reducing state
1336 calc_reducing_state();
1337}
1338
1339void HelmholtzEOSMixtureBackend::update(CoolProp::input_pairs input_pair, double value1, double value2) {
1340 if (get_debug_level() > 10) {
1341 std::cout << format("%s (%d): update called with (%d: (%s), %g, %g)", __FILE__, __LINE__, input_pair,
1342 get_input_pair_short_desc(input_pair).c_str(), value1, value2)
1343 << std::endl;
1344 }
1345
1346 CoolPropDbl ld_value1 = value1, ld_value2 = value2;
1347 pre_update(input_pair, ld_value1, ld_value2);
1348 value1 = ld_value1;
1349 value2 = ld_value2;
1350
1351 switch (input_pair) {
1352 case PT_INPUTS:
1353 _p = value1;
1354 _T = value2;
1355 FlashRoutines::PT_flash(*this);
1356 break;
1357 case DmolarT_INPUTS:
1358 _rhomolar = value1;
1359 _T = value2;
1360 FlashRoutines::DHSU_T_flash(*this, iDmolar);
1361 break;
1362 case SmolarT_INPUTS:
1363 _smolar = value1;
1364 _T = value2;
1365 FlashRoutines::DHSU_T_flash(*this, iSmolar);
1366 break;
1367 //case HmolarT_INPUTS:
1368 // _hmolar = value1; _T = value2; FlashRoutines::DHSU_T_flash(*this, iHmolar); break;
1369 //case TUmolar_INPUTS:
1370 // _T = value1; _umolar = value2; FlashRoutines::DHSU_T_flash(*this, iUmolar); break;
1371 case DmolarP_INPUTS:
1372 _rhomolar = value1;
1373 _p = value2;
1374 FlashRoutines::DP_flash(*this);
1375 break;
1376 case DmolarHmolar_INPUTS:
1377 _rhomolar = value1;
1378 _hmolar = value2;
1379 FlashRoutines::HSU_D_flash(*this, iHmolar);
1380 break;
1381 case DmolarSmolar_INPUTS:
1382 _rhomolar = value1;
1383 _smolar = value2;
1384 FlashRoutines::HSU_D_flash(*this, iSmolar);
1385 break;
1386 case DmolarUmolar_INPUTS:
1387 _rhomolar = value1;
1388 _umolar = value2;
1389 FlashRoutines::HSU_D_flash(*this, iUmolar);
1390 break;
1391 case HmolarP_INPUTS:
1392 _hmolar = value1;
1393 _p = value2;
1394 FlashRoutines::HSU_P_flash(*this, iHmolar);
1395 break;
1396 case PSmolar_INPUTS:
1397 _p = value1;
1398 _smolar = value2;
1399 FlashRoutines::HSU_P_flash(*this, iSmolar);
1400 break;
1401 case PUmolar_INPUTS:
1402 _p = value1;
1403 _umolar = value2;
1404 FlashRoutines::HSU_P_flash(*this, iUmolar);
1405 break;
1406 case HmolarSmolar_INPUTS:
1407 _hmolar = value1;
1408 _smolar = value2;
1409 FlashRoutines::HS_flash(*this);
1410 break;
1411 case QT_INPUTS:
1412 _Q = value1;
1413 _T = value2;
1414 if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1415 FlashRoutines::QT_flash(*this);
1416 break;
1417 case PQ_INPUTS:
1418 _p = value1;
1419 _Q = value2;
1420 if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1421 FlashRoutines::PQ_flash(*this);
1422 break;
1423 case QSmolar_INPUTS:
1424 _Q = value1;
1425 _smolar = value2;
1426 if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1427 FlashRoutines::QS_flash(*this);
1428 break;
1429 case HmolarQ_INPUTS:
1430 _hmolar = value1;
1431 _Q = value2;
1432 if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1433 FlashRoutines::HQ_flash(*this);
1434 break;
1435 case DmolarQ_INPUTS:
1436 _rhomolar = value1;
1437 _Q = value2;
1438 if ((_Q < 0) || (_Q > 1)) throw CoolProp::OutOfRangeError("Input vapor quality [Q] must be between 0 and 1");
1439 FlashRoutines::DQ_flash(*this);
1440 break;
1441 default:
1442 throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
1443 }
1444
1445 post_update();
1446}
1447const std::vector<CoolPropDbl> HelmholtzEOSMixtureBackend::calc_mass_fractions() {
1448 // mass fraction is mass_i/total_mass;
1449 CoolPropDbl mm = molar_mass();
1450 std::vector<CoolPropDbl>& mole_fractions = get_mole_fractions_ref();
1451 std::vector<CoolPropDbl> mass_fractions(mole_fractions.size());
1452 for (std::size_t i = 0; i < mole_fractions.size(); ++i) {
1453 double mmi = get_fluid_constant(i, imolar_mass);
1454 mass_fractions[i] = mmi * (mole_fractions[i]) / mm;
1455 }
1456 return mass_fractions;
1457}
1458
1459void HelmholtzEOSMixtureBackend::update_with_guesses(CoolProp::input_pairs input_pair, double value1, double value2,
1460 const GuessesStructure& guesses) {
1461 if (get_debug_level() > 10) {
1462 std::cout << format("%s (%d): update called with (%d: (%s), %g, %g)", __FILE__, __LINE__, input_pair,
1463 get_input_pair_short_desc(input_pair).c_str(), value1, value2)
1464 << std::endl;
1465 }
1466
1467 CoolPropDbl ld_value1 = value1, ld_value2 = value2;
1468 pre_update(input_pair, ld_value1, ld_value2);
1469 value1 = ld_value1;
1470 value2 = ld_value2;
1471
1472 switch (input_pair) {
1473 case PQ_INPUTS:
1474 _p = value1;
1475 _Q = value2;
1476 FlashRoutines::PQ_flash_with_guesses(*this, guesses);
1477 break;
1478 case QT_INPUTS:
1479 _Q = value1;
1480 _T = value2;
1481 FlashRoutines::QT_flash_with_guesses(*this, guesses);
1482 break;
1483 case PT_INPUTS:
1484 _p = value1;
1485 _T = value2;
1486 FlashRoutines::PT_flash_with_guesses(*this, guesses);
1487 break;
1488 default:
1489 throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
1490 }
1491 post_update();
1492}
1493
1494void HelmholtzEOSMixtureBackend::post_update(bool optional_checks) {
1495 // Check the values that must always be set
1496 //if (_p < 0){ throw ValueError("p is less than zero");}
1497 if (!ValidNumber(_p)) {
1498 throw ValueError("p is not a valid number");
1499 }
1500 //if (_T < 0){ throw ValueError("T is less than zero");}
1501 if (!ValidNumber(_T)) {
1502 throw ValueError("T is not a valid number");
1503 }
1504 if (_rhomolar < 0) {
1505 throw ValueError("rhomolar is less than zero");
1506 }
1507 if (!ValidNumber(_rhomolar)) {
1508 throw ValueError("rhomolar is not a valid number");
1509 }
1510
1511 if (optional_checks) {
1512 if (!ValidNumber(_Q)) {
1513 throw ValueError("Q is not a valid number");
1514 }
1515 if (_phase == iphase_unknown) {
1516 throw ValueError("_phase is unknown");
1517 }
1518 }
1519
1520 // Set the reduced variables
1521 _tau = _reducing.T / _T;
1522 _delta = _rhomolar / _reducing.rhomolar;
1523
1524 // Update the terms in the excess contribution
1525 if (!is_pure_or_pseudopure){
1526 residual_helmholtz->Excess.update(_tau, _delta);
1527 }
1528}
1529
1530CoolPropDbl HelmholtzEOSMixtureBackend::calc_Bvirial() {
1531 return 1 / rhomolar_reducing() * calc_alphar_deriv_nocache(0, 1, mole_fractions, _tau, 1e-12);
1532}
1533CoolPropDbl HelmholtzEOSMixtureBackend::calc_dBvirial_dT() {
1534 SimpleState red = get_reducing_state();
1535 CoolPropDbl dtau_dT = -red.T / pow(_T, 2);
1536 return 1 / red.rhomolar * calc_alphar_deriv_nocache(1, 1, mole_fractions, _tau, 1e-12) * dtau_dT;
1537}
1538CoolPropDbl HelmholtzEOSMixtureBackend::calc_Cvirial() {
1539 return 1 / pow(rhomolar_reducing(), 2) * calc_alphar_deriv_nocache(0, 2, mole_fractions, _tau, 1e-12);
1540}
1541CoolPropDbl HelmholtzEOSMixtureBackend::calc_dCvirial_dT() {
1542 SimpleState red = get_reducing_state();
1543 CoolPropDbl dtau_dT = -red.T / pow(_T, 2);
1544 return 1 / pow(red.rhomolar, 2) * calc_alphar_deriv_nocache(1, 2, mole_fractions, _tau, 1e-12) * dtau_dT;
1545}
1546void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int other, CoolPropDbl value, bool& saturation_called) {
1547 /*
1548 Determine the phase given p and one other state variable
1549 */
1550 saturation_called = false;
1551
1552 // Reference declaration to save indexing
1553 CoolPropFluid& component = components[0];
1554
1555 // Maximum saturation temperature - Equal to critical pressure for pure fluids
1556 CoolPropDbl psat_max = calc_pmax_sat();
1557
1558 // Check supercritical pressure
1559 if (_p > psat_max) {
1560 _Q = 1e9;
1561 switch (other) {
1562 case iT: {
1563 if (_T > _crit.T) {
1564 this->_phase = iphase_supercritical;
1565 return;
1566 } else {
1567 this->_phase = iphase_supercritical_liquid;
1568 return;
1569 }
1570 }
1571 case iDmolar: {
1572 if (_rhomolar < _crit.rhomolar) {
1573 this->_phase = iphase_supercritical_gas;
1574 return;
1575 } else {
1576 this->_phase = iphase_supercritical_liquid;
1577 return;
1578 }
1579 }
1580 case iSmolar: {
1581 if (_smolar.pt() > _crit.smolar) {
1582 this->_phase = iphase_supercritical_gas;
1583 return;
1584 } else {
1585 this->_phase = iphase_supercritical_liquid;
1586 return;
1587 }
1588 }
1589 case iHmolar: {
1590 if (_hmolar.pt() > _crit.hmolar) {
1591 this->_phase = iphase_supercritical_gas;
1592 return;
1593 } else {
1594 this->_phase = iphase_supercritical_liquid;
1595 return;
1596 }
1597 }
1598 case iUmolar: {
1599 if (_umolar.pt() > _crit.umolar) {
1600 this->_phase = iphase_supercritical_gas;
1601 return;
1602 } else {
1603 this->_phase = iphase_supercritical_liquid;
1604 return;
1605 }
1606 }
1607 default: {
1608 throw ValueError("supercritical pressure but other invalid for now");
1609 }
1610 }
1611 }
1612 // Check between triple point pressure and psat_max
1613 else if (_p >= components[0].EOS().ptriple * 0.9999 && _p <= psat_max) {
1614 // First try the ancillaries, use them to determine the state if you can
1615
1616 // Calculate dew and bubble temps from the ancillaries (everything needs them)
1617 _TLanc = components[0].ancillaries.pL.invert(_p);
1618 _TVanc = components[0].ancillaries.pV.invert(_p);
1619
1620 bool definitely_two_phase = false;
1621
1622 // Try using the ancillaries for P,H,S if they are there
1623 switch (other) {
1624 case iT: {
1625
1626 if (has_melting_line()) {
1627 double Tm = melting_line(iT, iP, _p);
1628 if (get_config_bool(DONT_CHECK_PROPERTY_LIMITS)) {
1629 _phase = iphase_liquid;
1630 } else {
1631 if (_T < Tm - 0.001) {
1632 throw ValueError(format("For now, we don't support T [%g K] below Tmelt(p) [%g K]", _T, Tm));
1633 }
1634 }
1635 } else {
1636 if (get_config_bool(DONT_CHECK_PROPERTY_LIMITS)) {
1637 _phase = iphase_liquid;
1638 } else {
1639 if (_T < Tmin() - 0.001) {
1640 throw ValueError(format("For now, we don't support T [%g K] below Tmin(saturation) [%g K]", _T, Tmin()));
1641 }
1642 }
1643 }
1644
1645 CoolPropDbl T_vap = 0.1 + static_cast<double>(_TVanc);
1646 CoolPropDbl T_liq = -0.1 + static_cast<double>(_TLanc);
1647
1648 if (value > T_vap) {
1649 this->_phase = iphase_gas;
1650 _Q = -1000;
1651 return;
1652 } else if (value < T_liq) {
1653 this->_phase = iphase_liquid;
1654 _Q = 1000;
1655 return;
1656 }
1657 break;
1658 }
1659 case iHmolar: {
1660 if (!component.ancillaries.hL.enabled()) {
1661 break;
1662 }
1663 // Ancillaries are h-h_anchor, so add back h_anchor
1664 CoolPropDbl h_liq = component.ancillaries.hL.evaluate(_TLanc) + component.EOS().hs_anchor.hmolar;
1665 CoolPropDbl h_liq_error_band = component.ancillaries.hL.get_max_abs_error();
1666 CoolPropDbl h_vap = h_liq + component.ancillaries.hLV.evaluate(_TLanc);
1667 CoolPropDbl h_vap_error_band = h_liq_error_band + component.ancillaries.hLV.get_max_abs_error();
1668
1669 // HelmholtzEOSMixtureBackend HEOS(components);
1670 // HEOS.update(QT_INPUTS, 1, _TLanc);
1671 // double h1 = HEOS.hmolar();
1672 // HEOS.update(QT_INPUTS, 0, _TLanc);
1673 // double h0 = HEOS.hmolar();
1674
1675 // Check if in range given the accuracy of the fit
1676 if (value > h_vap + h_vap_error_band) {
1677 this->_phase = iphase_gas;
1678 _Q = -1000;
1679 return;
1680 } else if (value < h_liq - h_liq_error_band) {
1681 this->_phase = iphase_liquid;
1682 _Q = 1000;
1683 return;
1684 } else if (value > h_liq + h_liq_error_band && value < h_vap - h_vap_error_band) {
1685 definitely_two_phase = true;
1686 }
1687 break;
1688 }
1689 case iSmolar: {
1690 if (!component.ancillaries.sL.enabled()) {
1691 break;
1692 }
1693 // Ancillaries are s-s_anchor, so add back s_anchor
1694 CoolPropDbl s_anchor = component.EOS().hs_anchor.smolar;
1695 CoolPropDbl s_liq = component.ancillaries.sL.evaluate(_TLanc) + s_anchor;
1696 CoolPropDbl s_liq_error_band = component.ancillaries.sL.get_max_abs_error();
1697 CoolPropDbl s_vap = s_liq + component.ancillaries.sLV.evaluate(_TVanc);
1698 CoolPropDbl s_vap_error_band = s_liq_error_band + component.ancillaries.sLV.get_max_abs_error();
1699
1700 // Check if in range given the accuracy of the fit
1701 if (value > s_vap + s_vap_error_band) {
1702 this->_phase = iphase_gas;
1703 _Q = -1000;
1704 return;
1705 } else if (value < s_liq - s_liq_error_band) {
1706 this->_phase = iphase_liquid;
1707 _Q = 1000;
1708 return;
1709 } else if (value > s_liq + s_liq_error_band && value < s_vap - s_vap_error_band) {
1710 definitely_two_phase = true;
1711 }
1712 break;
1713 }
1714 case iUmolar: {
1715 if (!component.ancillaries.hL.enabled()) {
1716 break;
1717 }
1718 // u = h-p/rho
1719
1720 // Ancillaries are h-h_anchor, so add back h_anchor
1721 CoolPropDbl h_liq = component.ancillaries.hL.evaluate(_TLanc) + component.EOS().hs_anchor.hmolar;
1722 CoolPropDbl h_liq_error_band = component.ancillaries.hL.get_max_abs_error();
1723 CoolPropDbl h_vap = h_liq + component.ancillaries.hLV.evaluate(_TLanc);
1724 CoolPropDbl h_vap_error_band = h_liq_error_band + component.ancillaries.hLV.get_max_abs_error();
1725 CoolPropDbl rho_vap = component.ancillaries.rhoV.evaluate(_TVanc);
1726 CoolPropDbl rho_liq = component.ancillaries.rhoL.evaluate(_TLanc);
1727 CoolPropDbl u_liq = h_liq - _p / rho_liq;
1728 CoolPropDbl u_vap = h_vap - _p / rho_vap;
1729 CoolPropDbl u_liq_error_band = 1.5 * h_liq_error_band; // Most of error is in enthalpy
1730 CoolPropDbl u_vap_error_band = 1.5 * h_vap_error_band; // Most of error is in enthalpy
1731
1732 // Check if in range given the accuracy of the fit
1733 if (value > u_vap + u_vap_error_band) {
1734 this->_phase = iphase_gas;
1735 _Q = -1000;
1736 return;
1737 } else if (value < u_liq - u_liq_error_band) {
1738 this->_phase = iphase_liquid;
1739 _Q = 1000;
1740 return;
1741 } else if (value > u_liq + u_liq_error_band && value < u_vap - u_vap_error_band) {
1742 definitely_two_phase = true;
1743 }
1744 break;
1745 }
1746 default: {
1747 }
1748 }
1749
1750 // Now either density is an input, or an ancillary for h,s,u is missing
1751 // Always calculate the densities using the ancillaries
1752 if (!definitely_two_phase) {
1753 _rhoVanc = component.ancillaries.rhoV.evaluate(_TVanc);
1754 _rhoLanc = component.ancillaries.rhoL.evaluate(_TLanc);
1755 CoolPropDbl rho_vap = 0.95 * static_cast<double>(_rhoVanc);
1756 CoolPropDbl rho_liq = 1.05 * static_cast<double>(_rhoLanc);
1757 switch (other) {
1758 case iDmolar: {
1759 if (value < rho_vap) {
1760 this->_phase = iphase_gas;
1761 return;
1762 } else if (value > rho_liq) {
1763 this->_phase = iphase_liquid;
1764 return;
1765 }
1766 break;
1767 }
1768 }
1769 }
1770
1771 if (!is_pure_or_pseudopure) {
1772 throw ValueError("possibly two-phase inputs not supported for mixtures for now");
1773 }
1774
1775 // Actually have to use saturation information sadly
1776 // For the given pressure, find the saturation state
1777 // Run the saturation routines to determine the saturation densities and pressures
1778 HelmholtzEOSMixtureBackend HEOS(components);
1779 HEOS._p = this->_p;
1780 HEOS._Q = 0; // ?? What is the best to do here? Doesn't matter for our purposes since pure fluid
1781 FlashRoutines::PQ_flash(HEOS);
1782
1783 // We called the saturation routines, so HEOS.SatL and HEOS.SatV are now updated
1784 // with the saturated liquid and vapor values, which can therefore be used in
1785 // the other solvers
1786 saturation_called = true;
1787
1788 CoolPropDbl Q;
1789
1790 if (other == iT) {
1791 if (value < HEOS.SatL->T() - 100 * DBL_EPSILON) {
1792 this->_phase = iphase_liquid;
1793 _Q = -1000;
1794 return;
1795 } else if (value > HEOS.SatV->T() + 100 * DBL_EPSILON) {
1796 this->_phase = iphase_gas;
1797 _Q = 1000;
1798 return;
1799 } else {
1800 this->_phase = iphase_twophase;
1801 }
1802 }
1803 switch (other) {
1804 case iDmolar:
1805 Q = (1 / value - 1 / HEOS.SatL->rhomolar()) / (1 / HEOS.SatV->rhomolar() - 1 / HEOS.SatL->rhomolar());
1806 break;
1807 case iSmolar:
1808 Q = (value - HEOS.SatL->smolar()) / (HEOS.SatV->smolar() - HEOS.SatL->smolar());
1809 break;
1810 case iHmolar:
1811 Q = (value - HEOS.SatL->hmolar()) / (HEOS.SatV->hmolar() - HEOS.SatL->hmolar());
1812 break;
1813 case iUmolar:
1814 Q = (value - HEOS.SatL->umolar()) / (HEOS.SatV->umolar() - HEOS.SatL->umolar());
1815 break;
1816 default:
1817 throw ValueError(format("bad input for other"));
1818 }
1819 // TODO: Check the speed penalty of these calls
1820 // Update the states
1821 if (this->SatL) this->SatL->update(DmolarT_INPUTS, HEOS.SatL->rhomolar(), HEOS.SatL->T());
1822 if (this->SatV) this->SatV->update(DmolarT_INPUTS, HEOS.SatV->rhomolar(), HEOS.SatV->T());
1823 // Update the two-Phase variables
1824 _rhoLmolar = HEOS.SatL->rhomolar();
1825 _rhoVmolar = HEOS.SatV->rhomolar();
1826
1827 //
1828 if (Q < -1e-9) {
1829 this->_phase = iphase_liquid;
1830 _Q = -1000;
1831 return;
1832 } else if (Q > 1 + 1e-9) {
1833 this->_phase = iphase_gas;
1834 _Q = 1000;
1835 return;
1836 } else {
1837 this->_phase = iphase_twophase;
1838 }
1839
1840 _Q = Q;
1841 // Load the outputs
1842 _T = _Q * HEOS.SatV->T() + (1 - _Q) * HEOS.SatL->T();
1843 _rhomolar = 1 / (_Q / HEOS.SatV->rhomolar() + (1 - _Q) / HEOS.SatL->rhomolar());
1844 return;
1845 } else if (_p < components[0].EOS().ptriple * 0.9999) {
1846 if (other == iT) {
1847 if (_T > std::max(Tmin(), Ttriple())) {
1848 _phase = iphase_gas;
1849 } else {
1850 if (get_config_bool(DONT_CHECK_PROPERTY_LIMITS)) {
1851 _phase = iphase_gas;
1852 } else {
1853 throw NotImplementedError(format("For now, we don't support p [%g Pa] below ptriple [%g Pa] when T [%g] is less than Tmin [%g]",
1854 _p, components[0].EOS().ptriple, _T, std::max(Tmin(), Ttriple())));
1855 }
1856 }
1857 } else {
1858 _phase = iphase_gas;
1859 }
1860 } else {
1861 throw ValueError(format("The pressure [%g Pa] cannot be used in p_phase_determination", _p));
1862 }
1863}
1864void HelmholtzEOSMixtureBackend::calc_ssat_max(void) {
1865 class Residual : public FuncWrapper1D
1866 {
1867 public:
1868 HelmholtzEOSMixtureBackend* HEOS;
1869 Residual(HelmholtzEOSMixtureBackend& HEOS) : HEOS(&HEOS){};
1870 double call(double T) {
1871 HEOS->update(QT_INPUTS, 1, T);
1872 // dTdp_along_sat
1873 double dTdp_along_sat =
1874 HEOS->T() * (1 / HEOS->SatV->rhomolar() - 1 / HEOS->SatL->rhomolar()) / (HEOS->SatV->hmolar() - HEOS->SatL->hmolar());
1875 // dsdT_along_sat;
1876 return HEOS->SatV->first_partial_deriv(iSmolar, iT, iP) + HEOS->SatV->first_partial_deriv(iSmolar, iP, iT) / dTdp_along_sat;
1877 }
1878 };
1879 if (!ssat_max.is_valid() && ssat_max.exists != SsatSimpleState::SSAT_MAX_DOESNT_EXIST) {
1880 shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS_copy(new CoolProp::HelmholtzEOSMixtureBackend(get_components()));
1881 Residual resid(*HEOS_copy);
1882 const CoolProp::SimpleState& tripleV = HEOS_copy->get_components()[0].triple_vapor;
1883 double v1 = resid.call(hsat_max.T);
1884 double v2 = resid.call(tripleV.T);
1885 // If there is a sign change, there is a maxima, otherwise there is no local maxima/minima
1886 if (v1 * v2 < 0) {
1887 Brent(resid, hsat_max.T, tripleV.T, DBL_EPSILON, 1e-8, 30);
1888 ssat_max.T = resid.HEOS->T();
1889 ssat_max.p = resid.HEOS->p();
1890 ssat_max.rhomolar = resid.HEOS->rhomolar();
1891 ssat_max.hmolar = resid.HEOS->hmolar();
1892 ssat_max.smolar = resid.HEOS->smolar();
1893 ssat_max.exists = SsatSimpleState::SSAT_MAX_DOES_EXIST;
1894 } else {
1895 ssat_max.exists = SsatSimpleState::SSAT_MAX_DOESNT_EXIST;
1896 }
1897 }
1898}
1899void HelmholtzEOSMixtureBackend::calc_hsat_max(void) {
1900 class Residualhmax : public FuncWrapper1D
1901 {
1902 public:
1903 HelmholtzEOSMixtureBackend* HEOS;
1904 Residualhmax(HelmholtzEOSMixtureBackend& HEOS) : HEOS(&HEOS){};
1905 double call(double T) {
1906 HEOS->update(QT_INPUTS, 1, T);
1907 // dTdp_along_sat
1908 double dTdp_along_sat =
1909 HEOS->T() * (1 / HEOS->SatV->rhomolar() - 1 / HEOS->SatL->rhomolar()) / (HEOS->SatV->hmolar() - HEOS->SatL->hmolar());
1910 // dhdT_along_sat;
1911 return HEOS->SatV->first_partial_deriv(iHmolar, iT, iP) + HEOS->SatV->first_partial_deriv(iHmolar, iP, iT) / dTdp_along_sat;
1912 }
1913 };
1914 if (!hsat_max.is_valid()) {
1915 shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS_copy(new CoolProp::HelmholtzEOSMixtureBackend(get_components()));
1916 Residualhmax residhmax(*HEOS_copy);
1917 Brent(residhmax, T_critical() - 0.1, HEOS_copy->Ttriple() + 1, DBL_EPSILON, 1e-8, 30);
1918 hsat_max.T = residhmax.HEOS->T();
1919 hsat_max.p = residhmax.HEOS->p();
1920 hsat_max.rhomolar = residhmax.HEOS->rhomolar();
1921 hsat_max.hmolar = residhmax.HEOS->hmolar();
1922 hsat_max.smolar = residhmax.HEOS->smolar();
1923 }
1924}
1925void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int other, CoolPropDbl value) {
1926 if (!ValidNumber(value)) {
1927 throw ValueError(format("value to T_phase_determination_pure_or_pseudopure is invalid"));
1928 };
1929
1930 // T is known, another input P, T, H, S, U is given (all molar)
1931 if (_T < _crit.T && _p > _crit.p) {
1932 // Only ever true if (other = iP); otherwise _p = -HUGE
1933 _phase = iphase_supercritical_liquid;
1934 } else if (std::abs(_T - _crit.T) < 10 * DBL_EPSILON) // Exactly at Tcrit
1935 {
1936 switch (other) {
1937 case iDmolar:
1938 if (std::abs(_rhomolar - _crit.rhomolar) < 10 * DBL_EPSILON) {
1939 _phase = iphase_critical_point;
1940 break;
1941 } else if (_rhomolar > _crit.rhomolar) {
1942 _phase = iphase_supercritical_liquid;
1943 break;
1944 } else {
1945 _phase = iphase_supercritical_gas;
1946 break;
1947 }
1948 case iP: {
1949 if (std::abs(_p - _crit.p) < 10 * DBL_EPSILON) {
1950 _phase = iphase_critical_point;
1951 break;
1952 } else if (_p > _crit.p) {
1953 _phase = iphase_supercritical_liquid;
1954 break;
1955 } else {
1956 _phase = iphase_supercritical_gas;
1957 break;
1958 }
1959 }
1960 default:
1961 throw ValueError(format("T=Tcrit; invalid input for other to T_phase_determination_pure_or_pseudopure"));
1962 }
1963 } else if (_T < _crit.T) // Gas, 2-Phase, Liquid, or Supercritical Liquid Region
1964 {
1965 // Start to think about the saturation stuff
1966 // First try to use the ancillary equations if you are far enough away
1967 // You know how accurate the ancillary equations are thanks to using CoolProp code to refit them
1968 switch (other) {
1969 case iP: {
1970 _pLanc = components[0].ancillaries.pL.evaluate(_T);
1971 _pVanc = components[0].ancillaries.pV.evaluate(_T);
1972 CoolPropDbl p_vap = 0.98 * static_cast<double>(_pVanc);
1973 CoolPropDbl p_liq = 1.02 * static_cast<double>(_pLanc);
1974
1975 if (value < p_vap) {
1976 this->_phase = iphase_gas;
1977 _Q = -1000;
1978 return;
1979 } else if (value > p_liq) {
1980 this->_phase = iphase_liquid;
1981 _Q = 1000;
1982 return;
1983 } else if (!is_pure()) // pseudo-pure
1984 {
1985 // For pseudo-pure fluids, the ancillary pressure curves are the official
1986 // arbiter of the phase
1987 if (value > static_cast<CoolPropDbl>(_pLanc)) {
1988 this->_phase = iphase_liquid;
1989 _Q = 1000;
1990 return;
1991 } else if (value < static_cast<CoolPropDbl>(_pVanc)) {
1992 this->_phase = iphase_gas;
1993 _Q = -1000;
1994 return;
1995 } else {
1996 throw ValueError("Two-phase inputs not supported for pseudo-pure for now");
1997 }
1998 }
1999 break;
2000 }
2001 default: {
2002 // Always calculate the densities using the ancillaries
2003 _rhoVanc = components[0].ancillaries.rhoV.evaluate(_T);
2004 _rhoLanc = components[0].ancillaries.rhoL.evaluate(_T);
2005 CoolPropDbl rho_vap = 0.95 * static_cast<double>(_rhoVanc);
2006 CoolPropDbl rho_liq = 1.05 * static_cast<double>(_rhoLanc);
2007 switch (other) {
2008 case iDmolar: {
2009 if (value < rho_vap) {
2010 this->_phase = iphase_gas;
2011 return;
2012 } else if (value > rho_liq) {
2013 this->_phase = iphase_liquid;
2014 return;
2015 } else {
2016 // Next we check the vapor quality based on the ancillary values
2017 double Qanc = (1 / value - 1 / static_cast<double>(_rhoLanc))
2018 / (1 / static_cast<double>(_rhoVanc) - 1 / static_cast<double>(_rhoLanc));
2019 // If the vapor quality is significantly inside the two-phase zone, stop, we are definitely two-phase
2020 if (value > 0.95 * rho_liq || value < 1.05 * rho_vap) {
2021 // Definitely single-phase
2022 _phase = iphase_liquid; // Needed for direct update call
2023 _Q = -1000; // Needed for direct update call
2024 update_DmolarT_direct(value, _T);
2025 CoolPropDbl pL = components[0].ancillaries.pL.evaluate(_T);
2026 if (Qanc < 0.01 && _p > pL * 1.05 && first_partial_deriv(iP, iDmolar, iT) > 0
2027 && second_partial_deriv(iP, iDmolar, iT, iDmolar, iT) > 0) {
2028 _phase = iphase_liquid;
2029 _Q = -1000;
2030 return;
2031 } else if (Qanc > 1.01) {
2032 break;
2033 } else {
2034 _phase = iphase_unknown;
2035 _p = _HUGE;
2036 }
2037 }
2038 }
2039 break;
2040 }
2041 default: {
2042 if (!this->SatL || !this->SatV) {
2043 throw ValueError(format("The saturation properties are needed in T_phase_determination_pure_or_pseudopure"));
2044 }
2045 // If it is not density, update the states
2046 SatV->update(DmolarT_INPUTS, rho_vap, _T);
2047 SatL->update(DmolarT_INPUTS, rho_liq, _T);
2048
2049 // First we check ancillaries
2050 switch (other) {
2051 case iSmolar: {
2052 if (value > SatV->calc_smolar()) {
2053 this->_phase = iphase_gas;
2054 return;
2055 }
2056 if (value < SatL->calc_smolar()) {
2057 this->_phase = iphase_liquid;
2058 return;
2059 }
2060 break;
2061 }
2062 case iHmolar: {
2063 if (value > SatV->calc_hmolar()) {
2064 this->_phase = iphase_gas;
2065 return;
2066 } else if (value < SatL->calc_hmolar()) {
2067 this->_phase = iphase_liquid;
2068 return;
2069 }
2070 break;
2071 }
2072 case iUmolar: {
2073 if (value > SatV->calc_umolar()) {
2074 this->_phase = iphase_gas;
2075 return;
2076 } else if (value < SatL->calc_umolar()) {
2077 this->_phase = iphase_liquid;
2078 return;
2079 }
2080 break;
2081 }
2082 default:
2083 throw ValueError(format("invalid input for other to T_phase_determination_pure_or_pseudopure"));
2084 }
2085 }
2086 }
2087 }
2088 }
2089
2090 // Actually have to use saturation information sadly
2091 // For the given temperature, find the saturation state
2092 // Run the saturation routines to determine the saturation densities and pressures
2093 HelmholtzEOSMixtureBackend HEOS(components);
2094 SaturationSolvers::saturation_T_pure_options options;
2095 SaturationSolvers::saturation_T_pure(HEOS, _T, options);
2096
2097 CoolPropDbl Q;
2098
2099 if (other == iP) {
2100 if (value > HEOS.SatL->p() * (1e-6 + 1)) {
2101 this->_phase = iphase_liquid;
2102 _Q = -1000;
2103 return;
2104 } else if (value < HEOS.SatV->p() * (1 - 1e-6)) {
2105 this->_phase = iphase_gas;
2106 _Q = 1000;
2107 return;
2108 } else {
2109 throw ValueError(
2110 format("Saturation pressure [%g Pa] corresponding to T [%g K] is within 1e-4 %% of given p [%Lg Pa]", HEOS.SatL->p(), _T, value));
2111 }
2112 }
2113
2114 switch (other) {
2115 case iDmolar:
2116 Q = (1 / value - 1 / HEOS.SatL->rhomolar()) / (1 / HEOS.SatV->rhomolar() - 1 / HEOS.SatL->rhomolar());
2117 break;
2118 case iSmolar:
2119 Q = (value - HEOS.SatL->smolar()) / (HEOS.SatV->smolar() - HEOS.SatL->smolar());
2120 break;
2121 case iHmolar:
2122 Q = (value - HEOS.SatL->hmolar()) / (HEOS.SatV->hmolar() - HEOS.SatL->hmolar());
2123 break;
2124 case iUmolar:
2125 Q = (value - HEOS.SatL->umolar()) / (HEOS.SatV->umolar() - HEOS.SatL->umolar());
2126 break;
2127 default:
2128 throw ValueError(format("bad input for other"));
2129 }
2130
2131 // Update the states
2132 if (this->SatL) this->SatL->update(DmolarT_INPUTS, HEOS.SatL->rhomolar(), HEOS.SatL->T());
2133 if (this->SatV) this->SatV->update(DmolarT_INPUTS, HEOS.SatV->rhomolar(), HEOS.SatV->T());
2134 // Update the two-Phase variables
2135 _rhoLmolar = HEOS.SatL->rhomolar();
2136 _rhoVmolar = HEOS.SatV->rhomolar();
2137
2138 if (Q < 0) {
2139 this->_phase = iphase_liquid;
2140 _Q = -1;
2141 return;
2142 } else if (Q > 1) {
2143 this->_phase = iphase_gas;
2144 _Q = 1;
2145 return;
2146 } else {
2147 this->_phase = iphase_twophase;
2148 }
2149 _Q = Q;
2150 // Load the outputs
2151 _p = _Q * HEOS.SatV->p() + (1 - _Q) * HEOS.SatL->p();
2152 _rhomolar = 1 / (_Q / HEOS.SatV->rhomolar() + (1 - _Q) / HEOS.SatL->rhomolar());
2153 return;
2154 } else if (_T > _crit.T && _T > components[0].EOS().Ttriple) // Supercritical or Supercritical Gas Region
2155 {
2156 _Q = 1e9;
2157 switch (other) {
2158 case iP: {
2159 if (_p > _crit.p) {
2160 this->_phase = iphase_supercritical;
2161 return;
2162 } else {
2163 this->_phase = iphase_supercritical_gas;
2164 return;
2165 }
2166 }
2167 case iDmolar: {
2168 if (_rhomolar > _crit.rhomolar) {
2169 this->_phase = iphase_supercritical_liquid;
2170 return;
2171 } else {
2172 this->_phase = iphase_supercritical_gas;
2173 return;
2174 }
2175 }
2176 case iSmolar: {
2177 if (_smolar.pt() > _crit.smolar) {
2178 this->_phase = iphase_supercritical_gas;
2179 return;
2180 } else {
2181 this->_phase = iphase_supercritical_liquid;
2182 return;
2183 }
2184 }
2185 case iHmolar: {
2186 if (_hmolar.pt() > _crit.hmolar) {
2187 this->_phase = iphase_supercritical_gas;
2188 return;
2189 } else {
2190 this->_phase = iphase_supercritical_liquid;
2191 return;
2192 }
2193 }
2194 case iUmolar: {
2195 if (_umolar.pt() > _crit.umolar) {
2196 this->_phase = iphase_supercritical_gas;
2197 return;
2198 } else {
2199 this->_phase = iphase_supercritical_liquid;
2200 return;
2201 }
2202 }
2203 default: {
2204 throw ValueError("supercritical temp but other invalid for now");
2205 }
2206 }
2207 } else {
2208 throw ValueError(format("For now, we don't support T [%g K] below Ttriple [%g K]", _T, components[0].EOS().Ttriple));
2209 }
2210}
2211void get_dT_drho(HelmholtzEOSMixtureBackend* HEOS, parameters index, CoolPropDbl& dT, CoolPropDbl& drho) {
2212 CoolPropDbl T = HEOS->T(), rho = HEOS->rhomolar(), rhor = HEOS->get_reducing_state().rhomolar, Tr = HEOS->get_reducing_state().T,
2213 dT_dtau = -pow(T, 2) / Tr, R = HEOS->gas_constant(), delta = rho / rhor, tau = Tr / T;
2214
2215 switch (index) {
2216 case iT:
2217 dT = 1;
2218 drho = 0;
2219 break;
2220 case iDmolar:
2221 dT = 0;
2222 drho = 1;
2223 break;
2224 case iDmass:
2225 dT = 0;
2226 drho = HEOS->molar_mass();
2227 break;
2228 case iP: {
2229 // dp/drho|T
2230 drho = R * T * (1 + 2 * delta * HEOS->dalphar_dDelta() + pow(delta, 2) * HEOS->d2alphar_dDelta2());
2231 // dp/dT|rho
2232 dT = rho * R * (1 + delta * HEOS->dalphar_dDelta() - tau * delta * HEOS->d2alphar_dDelta_dTau());
2233 break;
2234 }
2235 case iHmass:
2236 case iHmolar: {
2237 // dh/dT|rho
2238 dT = R
2239 * (-pow(tau, 2) * (HEOS->d2alpha0_dTau2() + HEOS->d2alphar_dTau2())
2240 + (1 + delta * HEOS->dalphar_dDelta() - tau * delta * HEOS->d2alphar_dDelta_dTau()));
2241 // dh/drhomolar|T
2242 drho =
2243 T * R / rho * (tau * delta * HEOS->d2alphar_dDelta_dTau() + delta * HEOS->dalphar_dDelta() + pow(delta, 2) * HEOS->d2alphar_dDelta2());
2244 if (index == iHmass) {
2245 // dhmolar/drhomolar|T * dhmass/dhmolar where dhmass/dhmolar = 1/mole mass
2246 drho /= HEOS->molar_mass();
2247 dT /= HEOS->molar_mass();
2248 }
2249 break;
2250 }
2251 case iSmass:
2252 case iSmolar: {
2253 // ds/dT|rho
2254 dT = R / T * (-pow(tau, 2) * (HEOS->d2alpha0_dTau2() + HEOS->d2alphar_dTau2()));
2255 // ds/drho|T
2256 drho = R / rho * (-(1 + delta * HEOS->dalphar_dDelta() - tau * delta * HEOS->d2alphar_dDelta_dTau()));
2257 if (index == iSmass) {
2258 // ds/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
2259 drho /= HEOS->molar_mass();
2260 dT /= HEOS->molar_mass();
2261 }
2262 break;
2263 }
2264 case iUmass:
2265 case iUmolar: {
2266 // du/dT|rho
2267 dT = R * (-pow(tau, 2) * (HEOS->d2alpha0_dTau2() + HEOS->d2alphar_dTau2()));
2268 // du/drho|T
2269 drho = HEOS->T() * R / rho * (tau * delta * HEOS->d2alphar_dDelta_dTau());
2270 if (index == iUmass) {
2271 // du/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
2272 drho /= HEOS->molar_mass();
2273 dT /= HEOS->molar_mass();
2274 }
2275 break;
2276 }
2277 case iTau:
2278 dT = 1 / dT_dtau;
2279 drho = 0;
2280 break;
2281 case iDelta:
2282 dT = 0;
2283 drho = 1 / rhor;
2284 break;
2285 default:
2286 throw ValueError(format("input to get_dT_drho[%s] is invalid", get_parameter_information(index, "short").c_str()));
2287 }
2288}
2289
2290CoolPropDbl HelmholtzEOSMixtureBackend::calc_pressure_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
2291 SimpleState reducing = calc_reducing_state_nocache(mole_fractions);
2292 CoolPropDbl delta = rhomolar / reducing.rhomolar;
2293 CoolPropDbl tau = reducing.T / T;
2294
2295 // Calculate derivative
2296 int nTau = 0, nDelta = 1;
2297 CoolPropDbl dalphar_dDelta = calc_alphar_deriv_nocache(nTau, nDelta, mole_fractions, tau, delta);
2298
2299 // Get pressure
2300 return rhomolar * gas_constant() * T * (1 + delta * dalphar_dDelta);
2301}
2302HelmholtzEOSBackend::StationaryPointReturnFlag HelmholtzEOSMixtureBackend::solver_dpdrho0_Tp(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomax,
2303 CoolPropDbl& light, CoolPropDbl& heavy) {
2304
2306 class dpdrho_resid : public FuncWrapper1DWithTwoDerivs
2307 {
2308 public:
2309 HelmholtzEOSMixtureBackend* HEOS;
2310 CoolPropDbl T, p, delta, rhor, tau, R_u;
2311
2312 dpdrho_resid(HelmholtzEOSMixtureBackend* HEOS, CoolPropDbl T, CoolPropDbl p)
2313 : HEOS(HEOS),
2314 T(T),
2315 p(p),
2316 delta(_HUGE),
2317 rhor(HEOS->get_reducing_state().rhomolar),
2318 tau(HEOS->get_reducing_state().T / T),
2319 R_u(HEOS->gas_constant()) {}
2320 double call(double rhomolar) {
2321 delta = rhomolar / rhor; // needed for derivative
2322 HEOS->update_DmolarT_direct(rhomolar, T);
2323 // dp/drho|T
2324 return R_u * T * (1 + 2 * delta * HEOS->dalphar_dDelta() + POW2(delta) * HEOS->d2alphar_dDelta2());
2325 };
2326 double deriv(double rhomolar) {
2327 // d2p/drho2|T
2328 return R_u * T / rhor * (2 * HEOS->dalphar_dDelta() + 4 * delta * HEOS->d2alphar_dDelta2() + POW2(delta) * HEOS->calc_d3alphar_dDelta3());
2329 };
2330 double second_deriv(double rhomolar) {
2331 // d3p/drho3|T
2332 return R_u * T / POW2(rhor)
2333 * (6 * HEOS->d2alphar_dDelta2() + 6 * delta * HEOS->d3alphar_dDelta3() + POW2(delta) * HEOS->calc_d4alphar_dDelta4());
2334 };
2335 };
2336 dpdrho_resid resid(this, T, p);
2337 light = -1;
2338 heavy = -1;
2339 try {
2340 light = Halley(resid, 1e-6, 1e-8, 100);
2341 double d2pdrho2__constT = resid.deriv(light);
2342 if (d2pdrho2__constT > 0) {
2343 // Not possible since curvature should be negative
2344 throw CoolProp::ValueError("curvature cannot be positive");
2345 }
2346 } catch (std::exception& e) {
2347 if (get_debug_level() > 5) {
2348 std::cout << e.what() << std::endl;
2349 };
2350 light = -1;
2351 }
2352
2353 if (light < 0) {
2354 try {
2355 // Now we are going to do something VERY slow - increase density until curvature is positive
2356 double rho = 1e-6;
2357 for (std::size_t counter = 0; counter <= 100; counter++) {
2358 resid.call(rho); // Updates the state
2359 double curvature = resid.deriv(rho);
2360 if (curvature > 0) {
2361 light = rho;
2362 break;
2363 }
2364 rho *= 2;
2365 }
2366 } catch (...) {
2367 }
2368 }
2369
2370 // First try a "normal" calculation of the stationary point on the liquid side
2371 for (double omega = 0.7; omega > 0; omega -= 0.2) {
2372 try {
2373 resid.options.add_number("omega", omega);
2374 heavy = Halley(resid, rhomax, 1e-8, 100);
2375 double d2pdrho2__constT = resid.deriv(heavy);
2376 if (d2pdrho2__constT < 0) {
2377 // Not possible since curvature should be positive
2378 throw CoolProp::ValueError("curvature cannot be negative");
2379 }
2380 break; // Jump out, we got a good solution
2381 } catch (std::exception& e) {
2382 if (get_debug_level() > 5) {
2383 std::cout << e.what() << std::endl;
2384 };
2385 heavy = -1;
2386 }
2387 }
2388
2389 if (heavy < 0) {
2390 try {
2391 // Now we are going to do something VERY slow - decrease density until curvature is negative or pressure is negative
2392 double rho = rhomax;
2393 for (std::size_t counter = 0; counter <= 100; counter++) {
2394 resid.call(rho); // Updates the state
2395 double curvature = resid.deriv(rho);
2396 if (curvature < 0 || this->p() < 0) {
2397 heavy = rho;
2398 break;
2399 }
2400 rho /= 1.1;
2401 }
2402 } catch (...) {
2403 }
2404 }
2405
2406 if (light > 0 && heavy > 0) {
2407 // Found two stationary points, done!
2408 return TWO_STATIONARY_POINTS_FOUND;
2409 }
2410 // If no solution is found for dpdrho|T=0 starting at high and low densities,
2411 // then try to do a bounded solver to see if you can find any solutions. If you
2412 // can't, p = f(rho) is probably monotonic (supercritical?), and the bounds are
2413 else if (light < 0 && heavy < 0) {
2414 double dpdrho_min = resid.call(1e-10);
2415 double dpdrho_max = resid.call(rhomax);
2416 if (dpdrho_max * dpdrho_min > 0) {
2417 return ZERO_STATIONARY_POINTS;
2418 } else {
2419 throw CoolProp::ValueError("zero stationary points -- does this make sense?");
2420 }
2421 } else {
2422 return ONE_STATIONARY_POINT_FOUND;
2423 }
2424}
2425// Define the residual to be driven to zero
2426class SolverTPResid : public FuncWrapper1DWithThreeDerivs
2427{
2428 public:
2429 HelmholtzEOSMixtureBackend* HEOS;
2430 CoolPropDbl T, p, delta, rhor, tau, R_u;
2431
2432 SolverTPResid(HelmholtzEOSMixtureBackend* HEOS, CoolPropDbl T, CoolPropDbl p)
2433 : HEOS(HEOS),
2434 T(T),
2435 p(p),
2436 delta(_HUGE),
2437 rhor(HEOS->get_reducing_state().rhomolar),
2438 tau(HEOS->get_reducing_state().T / T),
2439 R_u(HEOS->gas_constant()) {}
2440 double call(double rhomolar) {
2441 delta = rhomolar / rhor; // needed for derivative
2442 HEOS->update_DmolarT_direct(rhomolar, T);
2443 CoolPropDbl peos = HEOS->p();
2444 return (peos - p) / p;
2445 };
2446 double deriv(double rhomolar) {
2447 // dp/drho|T / pspecified
2448 return R_u * T * (1 + 2 * delta * HEOS->dalphar_dDelta() + POW2(delta) * HEOS->d2alphar_dDelta2()) / p;
2449 };
2450 double second_deriv(double rhomolar) {
2451 // d2p/drho2|T / pspecified
2452 return R_u * T / rhor * (2 * HEOS->dalphar_dDelta() + 4 * delta * HEOS->d2alphar_dDelta2() + POW2(delta) * HEOS->calc_d3alphar_dDelta3()) / p;
2453 };
2454 double third_deriv(double rhomolar) {
2455 // d3p/drho3|T / pspecified
2456 return R_u * T / POW2(rhor)
2457 * (6 * HEOS->d2alphar_dDelta2() + 6 * delta * HEOS->d3alphar_dDelta3() + POW2(delta) * HEOS->calc_d4alphar_dDelta4()) / p;
2458 };
2459};
2460CoolPropDbl HelmholtzEOSMixtureBackend::SRK_covolume() {
2461 double b = 0;
2462 for (std::size_t i = 0; i < mole_fractions.size(); ++i) {
2463 // Get the parameters for the cubic EOS
2464 CoolPropDbl Tc = get_fluid_constant(i, iT_critical), pc = get_fluid_constant(i, iP_critical);
2465 CoolPropDbl R = 8.3144598;
2466 b += mole_fractions[i] * 0.08664 * R * Tc / pc;
2467 }
2468 return b;
2469}
2470CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp_global(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomolar_max) {
2471 // Find the densities along the isotherm where dpdrho|T = 0 (if you can)
2472 CoolPropDbl light = -1, heavy = -1;
2473 StationaryPointReturnFlag retval = solver_dpdrho0_Tp(T, p, rhomolar_max, light, heavy);
2474
2475 // Define the solver class
2476 SolverTPResid resid(this, T, p);
2477
2478 if (retval == ZERO_STATIONARY_POINTS) {
2479 // It's monotonic (no stationary points found), so do the full bounded solver
2480 // for the density
2481 double rho = Brent(resid, 1e-10, rhomolar_max, DBL_EPSILON, 1e-8, 100);
2482 return rho;
2483 } else if (retval == TWO_STATIONARY_POINTS_FOUND) {
2484
2485 // Calculate the pressures at the min and max densities where dpdrho|T = 0
2486 double p_at_rhomin_stationary = calc_pressure_nocache(T, light);
2487 double p_at_rhomax_stationary = calc_pressure_nocache(T, heavy);
2488
2489 double rho_liq = -1, rho_vap = -1;
2490 if (p > p_at_rhomax_stationary) {
2491 int counter = 0;
2492 for (/* init above, for debugging */; counter <= 10; counter++) {
2493 // Bump up rhomax if needed to bound the given pressure
2494 double p_at_rhomax = calc_pressure_nocache(T, rhomolar_max);
2495 if (p_at_rhomax < p) {
2496 rhomolar_max *= 1.05;
2497 } else {
2498 break;
2499 }
2500 }
2501 // Look for liquid root starting at stationary point density
2502 rho_liq = Brent(resid, heavy, rhomolar_max, DBL_EPSILON, 1e-8, 100);
2503 }
2504
2505 if (p < p_at_rhomin_stationary) {
2506 // Look for vapor root starting at stationary point density
2507 rho_vap = Brent(resid, light, 1e-10, DBL_EPSILON, 1e-8, 100);
2508 }
2509
2510 if (rho_vap > 0 && rho_liq > 0) {
2511 // Both densities are the same
2512 if (std::abs(rho_vap - rho_liq) < 1e-10) {
2513 // return one of them
2514 return rho_vap;
2515 } else {
2516 // Two solutions found, keep the one with lower Gibbs energy
2517 double gibbsmolar_vap = calc_gibbsmolar_nocache(T, rho_vap);
2518 double gibbsmolar_liq = calc_gibbsmolar_nocache(T, rho_liq);
2519 if (gibbsmolar_liq < gibbsmolar_vap) {
2520 return rho_liq;
2521 } else {
2522 return rho_vap;
2523 }
2524 }
2525 } else if (rho_vap < 0 && rho_liq > 0) {
2526 // Liquid root found, return it
2527 return rho_liq;
2528 } else if (rho_vap > 0 && rho_liq < 0) {
2529 // Vapor root found, return it
2530 return rho_vap;
2531 } else {
2532 throw CoolProp::ValueError(format("No density solutions for T=%g,p=%g,z=%s", T, p, vec_to_string(static_cast<std::vector<double>>(mole_fractions), "%0.12g").c_str()));
2533 }
2534 } else {
2535 throw CoolProp::ValueError(
2536 format("One stationary point (not good) for T=%g,p=%g,z=%s", T, p, vec_to_string(static_cast<std::vector<double>>(mole_fractions), "%0.12g").c_str()));
2537 }
2538};
2539
2540CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomolar_guess) {
2541 phases phase;
2542
2543 SolverTPResid resid(this, T, p);
2544
2545 // Check if the phase is imposed
2546 if (imposed_phase_index != iphase_not_imposed)
2547 // Use the imposed phase index
2548 phase = imposed_phase_index;
2549 else
2550 // Use the phase index in the class
2551 phase = _phase;
2552
2553 if (rhomolar_guess < 0) // Not provided
2554 {
2555 // Calculate a guess value using SRK equation of state
2556 rhomolar_guess = solver_rho_Tp_SRK(T, p, phase);
2557
2558 // A gas-like phase, ideal gas might not be the perfect model, but probably good enough
2559 if (phase == iphase_gas || phase == iphase_supercritical_gas || phase == iphase_supercritical) {
2560 if (rhomolar_guess < 0 || !ValidNumber(rhomolar_guess)) // If the guess is bad, probably high temperature, use ideal gas
2561 {
2562 rhomolar_guess = p / (gas_constant() * T);
2563 }
2564 } else if (phase == iphase_liquid) {
2565 double rhomolar;
2566 if (is_pure_or_pseudopure) {
2567 // It's liquid at subcritical pressure, we can use ancillaries as guess value
2568 CoolPropDbl _rhoLancval = static_cast<CoolPropDbl>(components[0].ancillaries.rhoL.evaluate(T));
2569 try {
2570 // First we try with Halley's method starting at saturated liquid
2571 rhomolar = Halley(resid, _rhoLancval, 1e-8, 100);
2572 if (!ValidNumber(rhomolar) || first_partial_deriv(iP, iDmolar, iT) < 0
2573 || second_partial_deriv(iP, iDmolar, iT, iDmolar, iT) < 0) {
2574 throw ValueError("Liquid density is invalid");
2575 }
2576 } catch (std::exception&) {
2577 // Next we try with a Brent method bounded solver since the function is 1-1
2578 rhomolar = Brent(resid, _rhoLancval * 0.9, _rhoLancval * 1.3, DBL_EPSILON, 1e-8, 100);
2579 if (!ValidNumber(rhomolar)) {
2580 throw ValueError();
2581 }
2582 }
2583 } else {
2584 // Try with 4th order Householder method starting at a very high density
2585 rhomolar = Householder4(&resid, 3 * rhomolar_reducing(), 1e-8, 100);
2586 }
2587 return rhomolar;
2588 } else if (phase == iphase_supercritical_liquid) {
2589 CoolPropDbl rhoLancval = static_cast<CoolPropDbl>(components[0].ancillaries.rhoL.evaluate(T));
2590 // Next we try with a Brent method bounded solver since the function is 1-1
2591 double rhomolar = Brent(resid, rhoLancval * 0.99, rhomolar_critical() * 4, DBL_EPSILON, 1e-8, 100);
2592 if (!ValidNumber(rhomolar)) {
2593 throw ValueError();
2594 }
2595 return rhomolar;
2596 }
2597 }
2598
2599 try {
2600 double rhomolar = Householder4(resid, rhomolar_guess, 1e-8, 20);
2601 if (!ValidNumber(rhomolar) || rhomolar < 0) {
2602 throw ValueError();
2603 }
2604 if (phase == iphase_liquid) {
2605 double dpdrho = first_partial_deriv(iP, iDmolar, iT);
2606 double d2pdrho2 = second_partial_deriv(iP, iDmolar, iT, iDmolar, iT);
2607 if (dpdrho < 0 || d2pdrho2 < 0) {
2608 // Try again with a larger density in order to end up at the right solution
2609 rhomolar = Householder4(resid, 3 * rhomolar_reducing(), 1e-8, 100);
2610 return rhomolar;
2611 }
2612 } else if (phase == iphase_gas) {
2613 double dpdrho = first_partial_deriv(iP, iDmolar, iT);
2614 double d2pdrho2 = second_partial_deriv(iP, iDmolar, iT, iDmolar, iT);
2615 if (dpdrho < 0 || d2pdrho2 > 0) {
2616 // Try again with a tiny density in order to end up at the right solution
2617 rhomolar = Householder4(resid, 1e-6, 1e-8, 100);
2618 return rhomolar;
2619 }
2620 }
2621 return rhomolar;
2622 } catch (std::exception& e) {
2623 if (phase == iphase_supercritical || phase == iphase_supercritical_gas) {
2624 double rhomolar = Brent(resid, 1e-10, 3 * rhomolar_reducing(), DBL_EPSILON, 1e-8, 100);
2625 return rhomolar;
2626 } else if (is_pure_or_pseudopure && T > T_critical()) {
2627 try {
2628 double rhomolar = Brent(resid, 1e-10, 5 * rhomolar_reducing(), DBL_EPSILON, 1e-8, 100);
2629 return rhomolar;
2630
2631 } catch (...) {
2632 double rhomolar = Householder4(resid, 3 * rhomolar_reducing(), 1e-8, 100);
2633 return rhomolar;
2634 }
2635 }
2636 throw ValueError(format("solver_rho_Tp was unable to find a solution for T=%10Lg, p=%10Lg, with guess value %10Lg with error: %s", T, p,
2637 rhomolar_guess, e.what()));
2638 }
2639}
2640CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp_SRK(CoolPropDbl T, CoolPropDbl p, phases phase) {
2641 CoolPropDbl rhomolar, R_u = gas_constant(), a = 0, b = 0, k_ij = 0;
2642
2643 for (std::size_t i = 0; i < components.size(); ++i) {
2644 CoolPropDbl Tci = components[i].EOS().reduce.T, pci = components[i].EOS().reduce.p, acentric_i = components[i].EOS().acentric;
2645 CoolPropDbl m_i = 0.480 + 1.574 * acentric_i - 0.176 * pow(acentric_i, 2);
2646 CoolPropDbl b_i = 0.08664 * R_u * Tci / pci;
2647 b += mole_fractions[i] * b_i;
2648
2649 CoolPropDbl a_i = 0.42747 * pow(R_u * Tci, 2) / pci * pow(1 + m_i * (1 - sqrt(T / Tci)), 2);
2650
2651 for (std::size_t j = 0; j < components.size(); ++j) {
2652 CoolPropDbl Tcj = components[j].EOS().reduce.T, pcj = components[j].EOS().reduce.p, acentric_j = components[j].EOS().acentric;
2653 CoolPropDbl m_j = 0.480 + 1.574 * acentric_j - 0.176 * pow(acentric_j, 2);
2654
2655 CoolPropDbl a_j = 0.42747 * pow(R_u * Tcj, 2) / pcj * pow(1 + m_j * (1 - sqrt(T / Tcj)), 2);
2656
2657 k_ij = 0;
2658 //if (i == j){
2659 // k_ij = 0;
2660 //}
2661 //else{
2662 // k_ij = 0;
2663 //}
2664
2665 a += mole_fractions[i] * mole_fractions[j] * sqrt(a_i * a_j) * (1 - k_ij);
2666 }
2667 }
2668
2669 CoolPropDbl A = a * p / pow(R_u * T, 2);
2670 CoolPropDbl B = b * p / (R_u * T);
2671
2672 //Solve the cubic for solutions for Z = p/(rho*R*T)
2673 double Z0, Z1, Z2;
2674 int Nsolns;
2675 solve_cubic(1, -1, A - B - B * B, -A * B, Nsolns, Z0, Z1, Z2);
2676
2677 // Determine the guess value
2678 if (Nsolns == 1) {
2679 rhomolar = p / (Z0 * R_u * T);
2680 } else {
2681 CoolPropDbl rhomolar0 = p / (Z0 * R_u * T);
2682 CoolPropDbl rhomolar1 = p / (Z1 * R_u * T);
2683 CoolPropDbl rhomolar2 = p / (Z2 * R_u * T);
2684
2685 // Check if only one solution is positive, return the solution if that is the case
2686 if (rhomolar0 > 0 && rhomolar1 <= 0 && rhomolar2 <= 0) {
2687 return rhomolar0;
2688 }
2689 if (rhomolar0 <= 0 && rhomolar1 > 0 && rhomolar2 <= 0) {
2690 return rhomolar1;
2691 }
2692 if (rhomolar0 <= 0 && rhomolar1 <= 0 && rhomolar2 > 0) {
2693 return rhomolar2;
2694 }
2695
2696 switch (phase) {
2697 case iphase_liquid:
2698 case iphase_supercritical_liquid:
2699 rhomolar = max3(rhomolar0, rhomolar1, rhomolar2);
2700 break;
2701 case iphase_gas:
2702 case iphase_supercritical_gas:
2703 case iphase_supercritical:
2704 rhomolar = min3(rhomolar0, rhomolar1, rhomolar2);
2705 break;
2706 default:
2707 throw ValueError("Bad phase to solver_rho_Tp_SRK");
2708 };
2709 }
2710 return rhomolar;
2711}
2712
2713CoolPropDbl HelmholtzEOSMixtureBackend::calc_pressure(void) {
2714 // Calculate the reducing parameters
2715 _delta = _rhomolar / _reducing.rhomolar;
2716 _tau = _reducing.T / _T;
2717
2718 // Calculate derivative if needed
2719 CoolPropDbl dar_dDelta = dalphar_dDelta();
2720 CoolPropDbl R_u = gas_constant();
2721
2722 // Get pressure
2723 _p = _rhomolar * R_u * _T * (1 + _delta.pt() * dar_dDelta);
2724
2725 //std::cout << format("p: %13.12f %13.12f %10.9f %10.9f %10.9f %10.9f %g\n",_T,_rhomolar,_tau,_delta,mole_fractions[0],dar_dDelta,_p);
2726 //if (_p < 0){
2727 // throw ValueError("Pressure is less than zero");
2728 //}
2729
2730 return static_cast<CoolPropDbl>(_p);
2731}
2732CoolPropDbl HelmholtzEOSMixtureBackend::calc_hmolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
2733 // Calculate the reducing parameters
2734 CoolPropDbl delta = rhomolar / _reducing.rhomolar;
2735 CoolPropDbl tau = _reducing.T / T;
2736
2737 // Calculate derivatives if needed, or just use cached values
2738 // Calculate derivative if needed
2739 CoolPropDbl dar_dDelta = calc_alphar_deriv_nocache(0, 1, mole_fractions, tau, delta);
2740 CoolPropDbl dar_dTau = calc_alphar_deriv_nocache(1, 0, mole_fractions, tau, delta);
2741 CoolPropDbl da0_dTau = calc_alpha0_deriv_nocache(1, 0, mole_fractions, tau, delta, _reducing.T, _reducing.rhomolar);
2742 CoolPropDbl R_u = gas_constant();
2743
2744 // Get molar enthalpy
2745 return R_u * T * (1 + tau * (da0_dTau + dar_dTau) + delta * dar_dDelta);
2746}
2747CoolPropDbl HelmholtzEOSMixtureBackend::calc_hmolar(void) {
2748 if (get_debug_level() >= 50)
2749 std::cout << format("HelmholtzEOSMixtureBackend::calc_hmolar: 2phase: %d T: %g rhomomolar: %g", isTwoPhase(), _T, _rhomolar) << std::endl;
2750 if (isTwoPhase()) {
2751 if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for the two-phase properties"));
2752 if (std::abs(_Q) < DBL_EPSILON) {
2753 _hmolar = SatL->hmolar();
2754 } else if (std::abs(_Q - 1) < DBL_EPSILON) {
2755 _hmolar = SatV->hmolar();
2756 } else {
2757 _hmolar = _Q * SatV->hmolar() + (1 - _Q) * SatL->hmolar();
2758 }
2759 return static_cast<CoolPropDbl>(_hmolar);
2760 } else if (isHomogeneousPhase()) {
2761 // Calculate the reducing parameters
2762 _delta = _rhomolar / _reducing.rhomolar;
2763 _tau = _reducing.T / _T;
2764
2765 // Calculate derivatives if needed, or just use cached values
2766 CoolPropDbl da0_dTau = dalpha0_dTau();
2767 CoolPropDbl dar_dTau = dalphar_dTau();
2768 CoolPropDbl dar_dDelta = dalphar_dDelta();
2769 CoolPropDbl R_u = gas_constant();
2770
2771 // Get molar enthalpy
2772 _hmolar = R_u * _T * (1 + _tau.pt() * (da0_dTau + dar_dTau) + _delta.pt() * dar_dDelta);
2773
2774 return static_cast<CoolPropDbl>(_hmolar);
2775 } else {
2776 throw ValueError(format("phase is invalid in calc_hmolar"));
2777 }
2778}
2779CoolPropDbl HelmholtzEOSMixtureBackend::calc_smolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
2780 // Calculate the reducing parameters
2781 CoolPropDbl delta = rhomolar / _reducing.rhomolar;
2782 CoolPropDbl tau = _reducing.T / T;
2783
2784 // Calculate derivatives if needed, or just use cached values
2785 // Calculate derivative if needed
2786 CoolPropDbl dar_dTau = calc_alphar_deriv_nocache(1, 0, mole_fractions, tau, delta);
2787 CoolPropDbl ar = calc_alphar_deriv_nocache(0, 0, mole_fractions, tau, delta);
2788 CoolPropDbl da0_dTau = calc_alpha0_deriv_nocache(1, 0, mole_fractions, tau, delta, _reducing.T, _reducing.rhomolar);
2789 CoolPropDbl a0 = calc_alpha0_deriv_nocache(0, 0, mole_fractions, tau, delta, _reducing.T, _reducing.rhomolar);
2790 CoolPropDbl R_u = gas_constant();
2791
2792 // Get molar entropy
2793 return R_u * (tau * (da0_dTau + dar_dTau) - a0 - ar);
2794}
2795CoolPropDbl HelmholtzEOSMixtureBackend::calc_smolar(void) {
2796 if (isTwoPhase()) {
2797 if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for the two-phase properties"));
2798 if (std::abs(_Q) < DBL_EPSILON) {
2799 _smolar = SatL->smolar();
2800 } else if (std::abs(_Q - 1) < DBL_EPSILON) {
2801 _smolar = SatV->smolar();
2802 } else {
2803 _smolar = _Q * SatV->smolar() + (1 - _Q) * SatL->smolar();
2804 }
2805 return static_cast<CoolPropDbl>(_smolar);
2806 } else if (isHomogeneousPhase()) {
2807 // Calculate the reducing parameters
2808 _delta = _rhomolar / _reducing.rhomolar;
2809 _tau = _reducing.T / _T;
2810
2811 // Calculate derivatives if needed, or just use cached values
2812 CoolPropDbl da0_dTau = dalpha0_dTau();
2813 CoolPropDbl ar = alphar();
2814 CoolPropDbl a0 = alpha0();
2815 CoolPropDbl dar_dTau = dalphar_dTau();
2816 CoolPropDbl R_u = gas_constant();
2817
2818 // Get molar entropy
2819 _smolar = R_u * (_tau.pt() * (da0_dTau + dar_dTau) - a0 - ar);
2820
2821 return static_cast<CoolPropDbl>(_smolar);
2822 } else {
2823 throw ValueError(format("phase is invalid in calc_smolar"));
2824 }
2825}
2826CoolPropDbl HelmholtzEOSMixtureBackend::calc_umolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
2827 // Calculate the reducing parameters
2828 CoolPropDbl delta = rhomolar / _reducing.rhomolar;
2829 CoolPropDbl tau = _reducing.T / T;
2830
2831 // Calculate derivatives
2832 CoolPropDbl dar_dTau = calc_alphar_deriv_nocache(1, 0, mole_fractions, tau, delta);
2833 CoolPropDbl da0_dTau = calc_alpha0_deriv_nocache(1, 0, mole_fractions, tau, delta, _reducing.T, _reducing.rhomolar);
2834 CoolPropDbl R_u = gas_constant();
2835
2836 // Get molar internal energy
2837 return R_u * T * tau * (da0_dTau + dar_dTau);
2838}
2839CoolPropDbl HelmholtzEOSMixtureBackend::calc_umolar(void) {
2840 if (isTwoPhase()) {
2841 if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for the two-phase properties"));
2842 if (std::abs(_Q) < DBL_EPSILON) {
2843 _umolar = SatL->umolar();
2844 } else if (std::abs(_Q - 1) < DBL_EPSILON) {
2845 _umolar = SatV->umolar();
2846 } else {
2847 _umolar = _Q * SatV->umolar() + (1 - _Q) * SatL->umolar();
2848 }
2849 return static_cast<CoolPropDbl>(_umolar);
2850 } else if (isHomogeneousPhase()) {
2851 // Calculate the reducing parameters
2852 _delta = _rhomolar / _reducing.rhomolar;
2853 _tau = _reducing.T / _T;
2854
2855 // Calculate derivatives if needed, or just use cached values
2856 CoolPropDbl da0_dTau = dalpha0_dTau();
2857 CoolPropDbl dar_dTau = dalphar_dTau();
2858 CoolPropDbl R_u = gas_constant();
2859
2860 // Get molar internal energy
2861 _umolar = R_u * _T * _tau.pt() * (da0_dTau + dar_dTau);
2862
2863 return static_cast<CoolPropDbl>(_umolar);
2864 } else {
2865 throw ValueError(format("phase is invalid in calc_umolar"));
2866 }
2867}
2868CoolPropDbl HelmholtzEOSMixtureBackend::calc_cvmolar(void) {
2869 // Calculate the reducing parameters
2870 _delta = _rhomolar / _reducing.rhomolar;
2871 _tau = _reducing.T / _T;
2872
2873 // Calculate derivatives if needed, or just use cached values
2874 CoolPropDbl d2ar_dTau2 = d2alphar_dTau2();
2875 CoolPropDbl d2a0_dTau2 = d2alpha0_dTau2();
2876 CoolPropDbl R_u = gas_constant();
2877
2878 // Get cv
2879 _cvmolar = -R_u * pow(_tau.pt(), 2) * (d2ar_dTau2 + d2a0_dTau2);
2880
2881 return static_cast<double>(_cvmolar);
2882}
2883CoolPropDbl HelmholtzEOSMixtureBackend::calc_cpmolar(void) {
2884 // Calculate the reducing parameters
2885 _delta = _rhomolar / _reducing.rhomolar;
2886 _tau = _reducing.T / _T;
2887
2888 // Calculate derivatives if needed, or just use cached values
2889 CoolPropDbl d2a0_dTau2 = d2alpha0_dTau2();
2890 CoolPropDbl dar_dDelta = dalphar_dDelta();
2891 CoolPropDbl d2ar_dDelta2 = d2alphar_dDelta2();
2892 CoolPropDbl d2ar_dDelta_dTau = d2alphar_dDelta_dTau();
2893 CoolPropDbl d2ar_dTau2 = d2alphar_dTau2();
2894 CoolPropDbl R_u = gas_constant();
2895
2896 // Get cp
2897 _cpmolar = R_u
2898 * (-pow(_tau.pt(), 2) * (d2ar_dTau2 + d2a0_dTau2)
2899 + pow(1 + _delta.pt() * dar_dDelta - _delta.pt() * _tau.pt() * d2ar_dDelta_dTau, 2)
2900 / (1 + 2 * _delta.pt() * dar_dDelta + pow(_delta.pt(), 2) * d2ar_dDelta2));
2901
2902 return static_cast<double>(_cpmolar);
2903}
2904CoolPropDbl HelmholtzEOSMixtureBackend::calc_cpmolar_idealgas(void) {
2905 // Calculate the reducing parameters
2906 _delta = _rhomolar / _reducing.rhomolar;
2907 _tau = _reducing.T / _T;
2908
2909 // Calculate derivatives if needed, or just use cached values
2910 CoolPropDbl d2a0_dTau2 = d2alpha0_dTau2();
2911 CoolPropDbl R_u = gas_constant();
2912
2913 // Get cp of the ideal gas
2914 return R_u * (1 + (-pow(_tau.pt(), 2)) * d2a0_dTau2);
2915}
2916CoolPropDbl HelmholtzEOSMixtureBackend::calc_speed_sound(void) {
2917 if (isTwoPhase()) {
2918 if (std::abs(_Q) < DBL_EPSILON) {
2919 return SatL->speed_sound();
2920 } else if (std::abs(_Q - 1) < DBL_EPSILON) {
2921 return SatV->speed_sound();
2922 } else {
2923 throw ValueError(format("Speed of sound is not defined for two-phase states because it depends on the distribution of phases."));
2924 }
2925 } else if (isHomogeneousPhase()) {
2926 // Calculate the reducing parameters
2927 _delta = _rhomolar / _reducing.rhomolar;
2928 _tau = _reducing.T / _T;
2929
2930 // Calculate derivatives if needed, or just use cached values
2931 CoolPropDbl d2a0_dTau2 = d2alpha0_dTau2();
2932 CoolPropDbl dar_dDelta = dalphar_dDelta();
2933 CoolPropDbl d2ar_dDelta2 = d2alphar_dDelta2();
2934 CoolPropDbl d2ar_dDelta_dTau = d2alphar_dDelta_dTau();
2935 CoolPropDbl d2ar_dTau2 = d2alphar_dTau2();
2936 CoolPropDbl R_u = gas_constant();
2937 CoolPropDbl mm = molar_mass();
2938
2939 // Get speed of sound
2940 _speed_sound = sqrt(
2941 R_u * _T / mm
2942 * (1 + 2 * _delta.pt() * dar_dDelta + pow(_delta.pt(), 2) * d2ar_dDelta2
2943 - pow(1 + _delta.pt() * dar_dDelta - _delta.pt() * _tau.pt() * d2ar_dDelta_dTau, 2) / (pow(_tau.pt(), 2) * (d2ar_dTau2 + d2a0_dTau2))));
2944
2945 return static_cast<CoolPropDbl>(_speed_sound);
2946 } else {
2947 throw ValueError(format("phase is invalid in calc_speed_sound"));
2948 }
2949}
2950
2951CoolPropDbl HelmholtzEOSMixtureBackend::calc_gibbsmolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar) {
2952 // Calculate the reducing parameters
2953 CoolPropDbl delta = rhomolar / _reducing.rhomolar;
2954 CoolPropDbl tau = _reducing.T / T;
2955
2956 // Calculate derivatives
2957 CoolPropDbl dar_dDelta = calc_alphar_deriv_nocache(0, 1, mole_fractions, tau, delta);
2958 CoolPropDbl ar = calc_alphar_deriv_nocache(0, 0, mole_fractions, tau, delta);
2959 CoolPropDbl a0 = calc_alpha0_deriv_nocache(0, 0, mole_fractions, tau, delta, _reducing.T, _reducing.rhomolar);
2960
2961 // Get molar gibbs function
2962 return gas_constant() * T * (1 + a0 + ar + delta * dar_dDelta);
2963}
2964CoolPropDbl HelmholtzEOSMixtureBackend::calc_gibbsmolar(void) {
2965 if (isTwoPhase()) {
2966 if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for the two-phase properties"));
2967 _gibbsmolar = _Q * SatV->gibbsmolar() + (1 - _Q) * SatL->gibbsmolar();
2968 return static_cast<CoolPropDbl>(_gibbsmolar);
2969 } else if (isHomogeneousPhase()) {
2970 // Calculate the reducing parameters
2971 _delta = _rhomolar / _reducing.rhomolar;
2972 _tau = _reducing.T / _T;
2973
2974 // Calculate derivatives if needed, or just use cached values
2975 CoolPropDbl ar = alphar();
2976 CoolPropDbl a0 = alpha0();
2977 CoolPropDbl dar_dDelta = dalphar_dDelta();
2978 CoolPropDbl R_u = gas_constant();
2979
2980 // Get molar gibbs function
2981 _gibbsmolar = R_u * _T * (1 + a0 + ar + _delta.pt() * dar_dDelta);
2982
2983 return static_cast<CoolPropDbl>(_gibbsmolar);
2984 } else {
2985 throw ValueError(format("phase is invalid in calc_gibbsmolar"));
2986 }
2987}
2988void HelmholtzEOSMixtureBackend::calc_excess_properties(void) {
2989 _gibbsmolar_excess = this->gibbsmolar(), _smolar_excess = this->smolar(), _hmolar_excess = this->hmolar();
2990 _umolar_excess = this->umolar();
2991 _volumemolar_excess = 1 / this->rhomolar();
2992 for (std::size_t i = 0; i < components.size(); ++i) {
2993 transient_pure_state.reset(new HelmholtzEOSBackend(components[i].name));
2994 transient_pure_state->update(PT_INPUTS, p(), T());
2995 double x_i = mole_fractions[i];
2996 double R = gas_constant();
2997 _gibbsmolar_excess = static_cast<double>(_gibbsmolar_excess) - x_i * (transient_pure_state->gibbsmolar() + R * T() * log(x_i));
2998 _hmolar_excess = static_cast<double>(_hmolar_excess) - x_i * transient_pure_state->hmolar();
2999 _umolar_excess = static_cast<double>(_umolar_excess) - x_i * transient_pure_state->umolar();
3000 _smolar_excess = static_cast<double>(_smolar_excess) - x_i * (transient_pure_state->smolar() - R * log(x_i));
3001 _volumemolar_excess = static_cast<double>(_volumemolar_excess) - x_i / transient_pure_state->rhomolar();
3002 }
3003 _helmholtzmolar_excess = static_cast<double>(_umolar_excess) - _T * static_cast<double>(_smolar_excess);
3004}
3005
3006CoolPropDbl HelmholtzEOSMixtureBackend::calc_helmholtzmolar(void) {
3007 if (isTwoPhase()) {
3008 if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for the two-phase properties"));
3009 _helmholtzmolar = _Q * SatV->helmholtzmolar() + (1 - _Q) * SatL->helmholtzmolar();
3010 return static_cast<CoolPropDbl>(_helmholtzmolar);
3011 } else if (isHomogeneousPhase()) {
3012 // Calculate the reducing parameters
3013 _delta = _rhomolar / _reducing.rhomolar;
3014 _tau = _reducing.T / _T;
3015
3016 // Calculate derivatives if needed, or just use cached values
3017 CoolPropDbl ar = alphar();
3018 CoolPropDbl a0 = alpha0();
3019 CoolPropDbl R_u = gas_constant();
3020
3021 // Get molar Helmholtz energy
3022 _helmholtzmolar = R_u * _T * (a0 + ar);
3023
3024 return static_cast<CoolPropDbl>(_helmholtzmolar);
3025 } else {
3026 throw ValueError(format("phase is invalid in calc_helmholtzmolar"));
3027 }
3028}
3029CoolPropDbl HelmholtzEOSMixtureBackend::calc_fugacity_coefficient(std::size_t i) {
3030 x_N_dependency_flag xN_flag = XN_DEPENDENT;
3031 return exp(MixtureDerivatives::ln_fugacity_coefficient(*this, i, xN_flag));
3032}
3033CoolPropDbl HelmholtzEOSMixtureBackend::calc_fugacity(std::size_t i) {
3034 x_N_dependency_flag xN_flag = XN_DEPENDENT;
3035 return MixtureDerivatives::fugacity_i(*this, i, xN_flag);
3036}
3037CoolPropDbl HelmholtzEOSMixtureBackend::calc_chemical_potential(std::size_t i) {
3038 x_N_dependency_flag xN_flag = XN_DEPENDENT;
3039 double Tci = get_fluid_constant(i, iT_critical);
3040 double rhoci = get_fluid_constant(i, irhomolar_critical);
3041 double dnar_dni__const_T_V_nj = MixtureDerivatives::dnalphar_dni__constT_V_nj(*this, i, xN_flag);
3042 double dna0_dni__const_T_V_nj =
3043 components[i].EOS().alpha0.base(tau() * (Tci / T_reducing()), delta() / (rhoci / rhomolar_reducing())) + 1 + log(mole_fractions[i]);
3044 return gas_constant() * T() * (dna0_dni__const_T_V_nj + dnar_dni__const_T_V_nj);
3045}
3046CoolPropDbl HelmholtzEOSMixtureBackend::calc_phase_identification_parameter(void) {
3047 return 2
3048 - rhomolar()
3049 * (second_partial_deriv(iP, iDmolar, iT, iT, iDmolar) / first_partial_deriv(iP, iT, iDmolar)
3050 - second_partial_deriv(iP, iDmolar, iT, iDmolar, iT) / first_partial_deriv(iP, iDmolar, iT));
3051}
3052
3053SimpleState HelmholtzEOSMixtureBackend::calc_reducing_state_nocache(const std::vector<CoolPropDbl>& mole_fractions) {
3054 SimpleState reducing;
3055 if (is_pure_or_pseudopure) {
3056 reducing = components[0].EOS().reduce;
3057 } else {
3058 reducing.T = Reducing->Tr(mole_fractions);
3059 reducing.rhomolar = Reducing->rhormolar(mole_fractions);
3060 }
3061 return reducing;
3062}
3063void HelmholtzEOSMixtureBackend::calc_reducing_state(void) {
3064 if (get_mole_fractions_ref().empty()) {
3065 throw ValueError("Mole fractions must be set before calling calc_reducing_state");
3066 }
3068 _reducing = calc_reducing_state_nocache(mole_fractions);
3069 _crit = _reducing;
3070}
3071void HelmholtzEOSMixtureBackend::calc_all_alphar_deriv_cache(const std::vector<CoolPropDbl>& mole_fractions, const CoolPropDbl& tau,
3072 const CoolPropDbl& delta) {
3073 deriv_counter++;
3074 bool cache_values = true;
3075 HelmholtzDerivatives derivs = residual_helmholtz->all(*this, get_mole_fractions_ref(), tau, delta, cache_values);
3076 _alphar = derivs.alphar;
3077 _dalphar_dDelta = derivs.dalphar_ddelta;
3078 _dalphar_dTau = derivs.dalphar_dtau;
3079 _d2alphar_dDelta2 = derivs.d2alphar_ddelta2;
3080 _d2alphar_dDelta_dTau = derivs.d2alphar_ddelta_dtau;
3081 _d2alphar_dTau2 = derivs.d2alphar_dtau2;
3082 _d3alphar_dDelta3 = derivs.d3alphar_ddelta3;
3083 _d3alphar_dDelta2_dTau = derivs.d3alphar_ddelta2_dtau;
3084 _d3alphar_dDelta_dTau2 = derivs.d3alphar_ddelta_dtau2;
3085 _d3alphar_dTau3 = derivs.d3alphar_dtau3;
3086 _d4alphar_dDelta4 = derivs.d4alphar_ddelta4;
3087 _d4alphar_dDelta3_dTau = derivs.d4alphar_ddelta3_dtau;
3088 _d4alphar_dDelta2_dTau2 = derivs.d4alphar_ddelta2_dtau2;
3089 _d4alphar_dDelta_dTau3 = derivs.d4alphar_ddelta_dtau3;
3090 _d4alphar_dTau4 = derivs.d4alphar_dtau4;
3091}
3092
3093CoolPropDbl HelmholtzEOSMixtureBackend::calc_alphar_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl>& mole_fractions,
3094 const CoolPropDbl& tau, const CoolPropDbl& delta) {
3095 bool cache_values = false;
3096 HelmholtzDerivatives derivs = residual_helmholtz->all(*this, mole_fractions, tau, delta, cache_values);
3097 return derivs.get(nTau, nDelta);
3098}
3099CoolPropDbl HelmholtzEOSMixtureBackend::calc_alpha0_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl>& mole_fractions,
3100 const CoolPropDbl& tau, const CoolPropDbl& delta, const CoolPropDbl& Tr,
3101 const CoolPropDbl& rhor) {
3102 CoolPropDbl val;
3103 if (components.size() == 0) {
3104 throw ValueError("No alpha0 derivatives are available");
3105 }
3106 if (is_pure_or_pseudopure) {
3107 EquationOfState& E = components[0].EOS();
3108 // In the case of cubics, we need to use the shifted tau^*=Tc/T and delta^*=rho/rhoc
3109 // rather than tau=Tr/T and delta=rho/rhor
3110 // For multiparameter EOS, this changes nothing because Tc/Tr = 1 and rhoc/rhor = 1
3111 double Tc = get_fluid_constant(0, iT_reducing), rhomolarc = get_fluid_constant(0, irhomolar_reducing);
3112
3113 // Cache the reducing temperature in some terms that need it (GERG-2004 models)
3114 E.alpha0.set_Tred(Tc);
3115 double taustar = Tc / Tr * tau, deltastar = rhor / rhomolarc * delta;
3116 if (nTau == 0 && nDelta == 0) {
3117 val = E.base0(taustar, deltastar);
3118 } else if (nTau == 0 && nDelta == 1) {
3119 val = E.dalpha0_dDelta(taustar, deltastar);
3120 } else if (nTau == 1 && nDelta == 0) {
3121 val = E.dalpha0_dTau(taustar, deltastar);
3122 } else if (nTau == 0 && nDelta == 2) {
3123 val = E.d2alpha0_dDelta2(taustar, deltastar);
3124 } else if (nTau == 1 && nDelta == 1) {
3125 val = E.d2alpha0_dDelta_dTau(taustar, deltastar);
3126 } else if (nTau == 2 && nDelta == 0) {
3127 val = E.d2alpha0_dTau2(taustar, deltastar);
3128 } else if (nTau == 0 && nDelta == 3) {
3129 val = E.d3alpha0_dDelta3(taustar, deltastar);
3130 } else if (nTau == 1 && nDelta == 2) {
3131 val = E.d3alpha0_dDelta2_dTau(taustar, deltastar);
3132 } else if (nTau == 2 && nDelta == 1) {
3133 val = E.d3alpha0_dDelta_dTau2(taustar, deltastar);
3134 } else if (nTau == 3 && nDelta == 0) {
3135 val = E.d3alpha0_dTau3(taustar, deltastar);
3136 } else {
3137 throw ValueError();
3138 }
3139 val *= pow(rhor / rhomolarc, nDelta);
3140 val /= pow(Tr / Tc, nTau);
3141 if (!ValidNumber(val)) {
3142 //calc_alpha0_deriv_nocache(nTau,nDelta,mole_fractions,tau,delta,Tr,rhor);
3143 throw ValueError(format("calc_alpha0_deriv_nocache returned invalid number with inputs nTau: %d, nDelta: %d, tau: %Lg, delta: %Lg", nTau,
3144 nDelta, tau, delta));
3145 } else {
3146 return val;
3147 }
3148 } else {
3149 // See Table B5, GERG 2008 from Kunz Wagner, JCED, 2012
3150 std::size_t N = mole_fractions.size();
3151 CoolPropDbl summer = 0;
3152 CoolPropDbl tau_i, delta_i, rho_ci, T_ci;
3153 CoolPropDbl Rmix = gas_constant();
3154 for (unsigned int i = 0; i < N; ++i) {
3155
3156 rho_ci = get_fluid_constant(i, irhomolar_critical);
3157 T_ci = get_fluid_constant(i, iT_critical);
3158 CoolPropDbl Rcomponent = get_fluid_constant(i, igas_constant);
3159 tau_i = T_ci * tau / Tr;
3160 delta_i = delta * rhor / rho_ci;
3161 CoolPropDbl Rratio = Rcomponent / Rmix;
3162
3163 // Cache the reducing temperature in some terms that need it (GERG-2004 models)
3164 components[i].EOS().alpha0.set_Tred(Tr);
3165
3166 if (nTau == 0 && nDelta == 0) {
3167 double logxi = (std::abs(mole_fractions[i]) > DBL_EPSILON) ? log(mole_fractions[i]) : 0;
3168 summer += mole_fractions[i] * Rratio * (components[i].EOS().base0(tau_i, delta_i) + logxi);
3169 } else if (nTau == 0 && nDelta == 1) {
3170 summer += mole_fractions[i] * Rratio * rhor / rho_ci * components[i].EOS().dalpha0_dDelta(tau_i, delta_i);
3171 } else if (nTau == 1 && nDelta == 0) {
3172 summer += mole_fractions[i] * Rratio * T_ci / Tr * components[i].EOS().dalpha0_dTau(tau_i, delta_i);
3173 } else if (nTau == 0 && nDelta == 2) {
3174 summer += mole_fractions[i] * Rratio * pow(rhor / rho_ci, 2) * components[i].EOS().d2alpha0_dDelta2(tau_i, delta_i);
3175 } else if (nTau == 1 && nDelta == 1) {
3176 summer += mole_fractions[i] * Rratio * rhor / rho_ci * T_ci / Tr * components[i].EOS().d2alpha0_dDelta_dTau(tau_i, delta_i);
3177 } else if (nTau == 2 && nDelta == 0) {
3178 summer += mole_fractions[i] * Rratio * pow(T_ci / Tr, 2) * components[i].EOS().d2alpha0_dTau2(tau_i, delta_i);
3179 } else {
3180 throw ValueError();
3181 }
3182 }
3183 return summer;
3184 }
3185}
3186CoolPropDbl HelmholtzEOSMixtureBackend::calc_alphar(void) {
3187 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3188 return static_cast<CoolPropDbl>(_alphar);
3189}
3190CoolPropDbl HelmholtzEOSMixtureBackend::calc_dalphar_dDelta(void) {
3191 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3192 return static_cast<CoolPropDbl>(_dalphar_dDelta);
3193}
3194CoolPropDbl HelmholtzEOSMixtureBackend::calc_dalphar_dTau(void) {
3195 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3196 return static_cast<CoolPropDbl>(_dalphar_dTau);
3197}
3198CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alphar_dTau2(void) {
3199 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3200 return static_cast<CoolPropDbl>(_d2alphar_dTau2);
3201}
3202CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alphar_dDelta_dTau(void) {
3203 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3204 return static_cast<CoolPropDbl>(_d2alphar_dDelta_dTau);
3205}
3206CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alphar_dDelta2(void) {
3207 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3208 return static_cast<CoolPropDbl>(_d2alphar_dDelta2);
3209}
3210CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alphar_dDelta3(void) {
3211 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3212 return static_cast<CoolPropDbl>(_d3alphar_dDelta3);
3213}
3214CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alphar_dDelta2_dTau(void) {
3215 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3216 return static_cast<CoolPropDbl>(_d3alphar_dDelta2_dTau);
3217}
3218CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alphar_dDelta_dTau2(void) {
3219 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3220 return static_cast<CoolPropDbl>(_d3alphar_dDelta_dTau2);
3221}
3222CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alphar_dTau3(void) {
3223 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3224 return static_cast<CoolPropDbl>(_d3alphar_dTau3);
3225}
3226
3227CoolPropDbl HelmholtzEOSMixtureBackend::calc_d4alphar_dDelta4(void) {
3228 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3229 return static_cast<CoolPropDbl>(_d4alphar_dDelta4);
3230}
3231CoolPropDbl HelmholtzEOSMixtureBackend::calc_d4alphar_dDelta3_dTau(void) {
3232 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3233 return static_cast<CoolPropDbl>(_d4alphar_dDelta3_dTau);
3234}
3235CoolPropDbl HelmholtzEOSMixtureBackend::calc_d4alphar_dDelta2_dTau2(void) {
3236 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3237 return static_cast<CoolPropDbl>(_d4alphar_dDelta2_dTau2);
3238}
3239CoolPropDbl HelmholtzEOSMixtureBackend::calc_d4alphar_dDelta_dTau3(void) {
3240 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3241 return static_cast<CoolPropDbl>(_d4alphar_dDelta_dTau3);
3242}
3243CoolPropDbl HelmholtzEOSMixtureBackend::calc_d4alphar_dTau4(void) {
3244 calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta);
3245 return static_cast<CoolPropDbl>(_d4alphar_dTau4);
3246}
3247
3248CoolPropDbl HelmholtzEOSMixtureBackend::calc_alpha0(void) {
3249 const int nTau = 0, nDelta = 0;
3250 return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3251}
3252CoolPropDbl HelmholtzEOSMixtureBackend::calc_dalpha0_dDelta(void) {
3253 const int nTau = 0, nDelta = 1;
3254 return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3255}
3256CoolPropDbl HelmholtzEOSMixtureBackend::calc_dalpha0_dTau(void) {
3257 const int nTau = 1, nDelta = 0;
3258 return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3259}
3260CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alpha0_dDelta2(void) {
3261 const int nTau = 0, nDelta = 2;
3262 return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3263}
3264CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alpha0_dDelta_dTau(void) {
3265 const int nTau = 1, nDelta = 1;
3266 return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3267}
3268CoolPropDbl HelmholtzEOSMixtureBackend::calc_d2alpha0_dTau2(void) {
3269 const int nTau = 2, nDelta = 0;
3270 return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3271}
3272CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alpha0_dDelta3(void) {
3273 const int nTau = 0, nDelta = 3;
3274 return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3275}
3276CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alpha0_dDelta2_dTau(void) {
3277 const int nTau = 1, nDelta = 2;
3278 return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3279}
3280CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alpha0_dDelta_dTau2(void) {
3281 const int nTau = 2, nDelta = 1;
3282 return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3283}
3284CoolPropDbl HelmholtzEOSMixtureBackend::calc_d3alpha0_dTau3(void) {
3285 const int nTau = 3, nDelta = 0;
3286 return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
3287}
3288CoolPropDbl HelmholtzEOSMixtureBackend::calc_first_saturation_deriv(parameters Of1, parameters Wrt1, HelmholtzEOSMixtureBackend& SatL,
3289 HelmholtzEOSMixtureBackend& SatV) {
3290 // Derivative of temperature w.r.t. pressure ALONG the saturation curve
3291 CoolPropDbl dTdP_sat = T() * (1 / SatV.rhomolar() - 1 / SatL.rhomolar()) / (SatV.hmolar() - SatL.hmolar());
3292
3293 // "Trivial" inputs
3294 if (Of1 == iT && Wrt1 == iP) {
3295 return dTdP_sat;
3296 } else if (Of1 == iP && Wrt1 == iT) {
3297 return 1 / dTdP_sat;
3298 }
3299 // Derivative taken with respect to T
3300 else if (Wrt1 == iT) {
3301 return first_partial_deriv(Of1, iT, iP) + first_partial_deriv(Of1, iP, iT) / dTdP_sat;
3302 }
3303 // Derivative taken with respect to p
3304 else if (Wrt1 == iP) {
3305 return first_partial_deriv(Of1, iP, iT) + first_partial_deriv(Of1, iT, iP) * dTdP_sat;
3306 } else {
3307 throw ValueError(
3308 format("Not possible to take first saturation derivative with respect to %s", get_parameter_information(Wrt1, "short").c_str()));
3309 }
3310}
3311CoolPropDbl HelmholtzEOSMixtureBackend::calc_first_saturation_deriv(parameters Of1, parameters Wrt1) {
3312 if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for calc_first_saturation_deriv"));
3313
3314 // Derivative of temperature w.r.t. pressure ALONG the saturation curve
3315 CoolPropDbl dTdP_sat = T() * (1 / SatV->rhomolar() - 1 / SatL->rhomolar()) / (SatV->hmolar() - SatL->hmolar());
3316
3317 // "Trivial" inputs
3318 if (Of1 == iT && Wrt1 == iP) {
3319 return dTdP_sat;
3320 } else if (Of1 == iP && Wrt1 == iT) {
3321 return 1 / dTdP_sat;
3322 }
3323 // Derivative taken with respect to T
3324 else if (Wrt1 == iT) {
3325 return first_partial_deriv(Of1, iT, iP) + first_partial_deriv(Of1, iP, iT) / dTdP_sat;
3326 }
3327 // Derivative taken with respect to p
3328 else if (Wrt1 == iP) {
3329 return first_partial_deriv(Of1, iP, iT) + first_partial_deriv(Of1, iT, iP) * dTdP_sat;
3330 } else {
3331 throw ValueError(
3332 format("Not possible to take first saturation derivative with respect to %s", get_parameter_information(Wrt1, "short").c_str()));
3333 }
3334}
3335CoolPropDbl HelmholtzEOSMixtureBackend::calc_second_saturation_deriv(parameters Of1, parameters Wrt1, parameters Wrt2) {
3336 if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for calc_second_saturation_deriv"));
3337 if (Wrt1 == iP && Wrt2 == iP) {
3338 CoolPropDbl dydT_constp = this->first_partial_deriv(Of1, iT, iP);
3339 CoolPropDbl d2ydTdp = this->second_partial_deriv(Of1, iT, iP, iP, iT);
3340 CoolPropDbl d2ydp2_constT = this->second_partial_deriv(Of1, iP, iT, iP, iT);
3341 CoolPropDbl d2ydT2_constp = this->second_partial_deriv(Of1, iT, iP, iT, iP);
3342
3343 CoolPropDbl dTdp_along_sat = calc_first_saturation_deriv(iT, iP);
3344 CoolPropDbl dvdrhoL = -1 / POW2(SatL->rhomolar());
3345 CoolPropDbl dvdrhoV = -1 / POW2(SatV->rhomolar());
3346 CoolPropDbl DELTAv = 1 / SatV->rhomolar() - 1 / SatL->rhomolar();
3347 CoolPropDbl dDELTAv_dT_constp = dvdrhoV * SatV->first_partial_deriv(iDmolar, iT, iP) - dvdrhoL * SatL->first_partial_deriv(iDmolar, iT, iP);
3348 CoolPropDbl dDELTAv_dp_constT = dvdrhoV * SatV->first_partial_deriv(iDmolar, iP, iT) - dvdrhoL * SatL->first_partial_deriv(iDmolar, iP, iT);
3349 CoolPropDbl DELTAh = SatV->hmolar() - SatL->hmolar();
3350 CoolPropDbl dDELTAh_dT_constp = SatV->first_partial_deriv(iHmolar, iT, iP) - SatL->first_partial_deriv(iHmolar, iT, iP);
3351 CoolPropDbl dDELTAh_dp_constT = SatV->first_partial_deriv(iHmolar, iP, iT) - SatL->first_partial_deriv(iHmolar, iP, iT);
3352 CoolPropDbl ddT_dTdp_along_sat_constp = (DELTAh * (_T * dDELTAv_dT_constp + DELTAv) - _T * DELTAv * dDELTAh_dT_constp) / POW2(DELTAh);
3353 CoolPropDbl ddp_dTdp_along_sat_constT = (DELTAh * (_T * dDELTAv_dp_constT) - _T * DELTAv * dDELTAh_dp_constT) / POW2(DELTAh);
3354
3355 double ddp_dydpsigma = d2ydp2_constT + dydT_constp * ddp_dTdp_along_sat_constT + d2ydTdp * dTdp_along_sat;
3356 double ddT_dydpsigma = d2ydTdp + dydT_constp * ddT_dTdp_along_sat_constp + d2ydT2_constp * dTdp_along_sat;
3357 return ddp_dydpsigma + ddT_dydpsigma * dTdp_along_sat;
3358 } else {
3359 throw ValueError(format("Currently, only possible to take second saturation derivative w.r.t. P (both times)"));
3360 }
3361}
3362
3363CoolPropDbl HelmholtzEOSMixtureBackend::calc_second_two_phase_deriv(parameters Of, parameters Wrt1, parameters Constant1, parameters Wrt2,
3364 parameters Constant2) {
3365 if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for calc_second_two_phase_deriv"));
3366
3367 if (Of == iDmolar
3368 && ((Wrt1 == iHmolar && Constant1 == iP && Wrt2 == iP && Constant2 == iHmolar)
3369 || (Wrt2 == iHmolar && Constant2 == iP && Wrt1 == iP && Constant1 == iHmolar))) {
3370 parameters h_key = iHmolar, rho_key = iDmolar, p_key = iP;
3371 // taking the derivative of (drho/dv)*(dv/dh|p) with respect to p with h constant
3372 CoolPropDbl dv_dh_constp = calc_first_two_phase_deriv(rho_key, h_key, p_key) / (-POW2(rhomolar()));
3373 CoolPropDbl drhomolar_dp__consth = first_two_phase_deriv(rho_key, p_key, h_key);
3374
3375 // Calculate the derivative of dvdh|p with respect to p at constant h
3376 CoolPropDbl dhL_dp_sat = SatL->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
3377 CoolPropDbl dhV_dp_sat = SatV->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
3378 CoolPropDbl drhoL_dp_sat = SatL->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
3379 CoolPropDbl drhoV_dp_sat = SatV->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
3380 CoolPropDbl numerator = 1 / SatV->keyed_output(rho_key) - 1 / SatL->keyed_output(rho_key);
3381 CoolPropDbl denominator = SatV->keyed_output(h_key) - SatL->keyed_output(h_key);
3382 CoolPropDbl dnumerator = -1 / POW2(SatV->keyed_output(rho_key)) * drhoV_dp_sat + 1 / POW2(SatL->keyed_output(rho_key)) * drhoL_dp_sat;
3383 CoolPropDbl ddenominator = dhV_dp_sat - dhL_dp_sat;
3384 CoolPropDbl d_dvdh_dp__consth = (denominator * dnumerator - numerator * ddenominator) / POW2(denominator);
3385 return -POW2(rhomolar()) * d_dvdh_dp__consth + dv_dh_constp * (-2 * rhomolar()) * drhomolar_dp__consth;
3386 } else if (Of == iDmass
3387 && ((Wrt1 == iHmass && Constant1 == iP && Wrt2 == iP && Constant2 == iHmass)
3388 || (Wrt2 == iHmass && Constant2 == iP && Wrt1 == iP && Constant1 == iHmass))) {
3389 parameters h_key = iHmass, rho_key = iDmass, p_key = iP;
3390 CoolPropDbl rho = keyed_output(rho_key);
3391 // taking the derivative of (drho/dv)*(dv/dh|p) with respect to p with h constant
3392 CoolPropDbl dv_dh_constp = calc_first_two_phase_deriv(rho_key, h_key, p_key) / (-POW2(rho));
3393 CoolPropDbl drho_dp__consth = first_two_phase_deriv(rho_key, p_key, h_key);
3394
3395 // Calculate the derivative of dvdh|p with respect to p at constant h
3396 CoolPropDbl dhL_dp_sat = SatL->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
3397 CoolPropDbl dhV_dp_sat = SatV->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
3398 CoolPropDbl drhoL_dp_sat = SatL->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
3399 CoolPropDbl drhoV_dp_sat = SatV->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
3400 CoolPropDbl numerator = 1 / SatV->keyed_output(rho_key) - 1 / SatL->keyed_output(rho_key);
3401 CoolPropDbl denominator = SatV->keyed_output(h_key) - SatL->keyed_output(h_key);
3402 CoolPropDbl dnumerator = -1 / POW2(SatV->keyed_output(rho_key)) * drhoV_dp_sat + 1 / POW2(SatL->keyed_output(rho_key)) * drhoL_dp_sat;
3403 CoolPropDbl ddenominator = dhV_dp_sat - dhL_dp_sat;
3404 CoolPropDbl d_dvdh_dp__consth = (denominator * dnumerator - numerator * ddenominator) / POW2(denominator);
3405 return -POW2(rho) * d_dvdh_dp__consth + dv_dh_constp * (-2 * rho) * drho_dp__consth;
3406 } else {
3407 throw ValueError();
3408 }
3409}
3410CoolPropDbl HelmholtzEOSMixtureBackend::calc_first_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant) {
3411 if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for calc_first_two_phase_deriv"));
3412 if (Of == iDmolar && Wrt == iHmolar && Constant == iP) {
3413 return -POW2(rhomolar()) * (1 / SatV->rhomolar() - 1 / SatL->rhomolar()) / (SatV->hmolar() - SatL->hmolar());
3414 } else if (Of == iDmass && Wrt == iHmass && Constant == iP) {
3415 return -POW2(rhomass()) * (1 / SatV->rhomass() - 1 / SatL->rhomass()) / (SatV->hmass() - SatL->hmass());
3416 } else if (Of == iDmolar && Wrt == iP && Constant == iHmolar) {
3417 // v = 1/rho; drhodv = -rho^2; dvdrho = -1/rho^2
3418 CoolPropDbl dvdrhoL = -1 / POW2(SatL->rhomolar());
3419 CoolPropDbl dvdrhoV = -1 / POW2(SatV->rhomolar());
3420 CoolPropDbl dvL_dp = dvdrhoL * SatL->calc_first_saturation_deriv(iDmolar, iP, *SatL, *SatV);
3421 CoolPropDbl dvV_dp = dvdrhoV * SatV->calc_first_saturation_deriv(iDmolar, iP, *SatL, *SatV);
3422 CoolPropDbl dhL_dp = SatL->calc_first_saturation_deriv(iHmolar, iP, *SatL, *SatV);
3423 CoolPropDbl dhV_dp = SatV->calc_first_saturation_deriv(iHmolar, iP, *SatL, *SatV);
3424 CoolPropDbl dxdp_h = (Q() * dhV_dp + (1 - Q()) * dhL_dp) / (SatL->hmolar() - SatV->hmolar());
3425 CoolPropDbl dvdp_h = dvL_dp + dxdp_h * (1 / SatV->rhomolar() - 1 / SatL->rhomolar()) + Q() * (dvV_dp - dvL_dp);
3426 return -POW2(rhomolar()) * dvdp_h;
3427 } else if (Of == iDmass && Wrt == iP && Constant == iHmass) {
3428 // v = 1/rho; drhodv = -rho^2; dvdrho = -1/rho^2
3429 CoolPropDbl dvdrhoL = -1 / POW2(SatL->rhomass());
3430 CoolPropDbl dvdrhoV = -1 / POW2(SatV->rhomass());
3431 CoolPropDbl dvL_dp = dvdrhoL * SatL->calc_first_saturation_deriv(iDmass, iP, *SatL, *SatV);
3432 CoolPropDbl dvV_dp = dvdrhoV * SatV->calc_first_saturation_deriv(iDmass, iP, *SatL, *SatV);
3433 CoolPropDbl dhL_dp = SatL->calc_first_saturation_deriv(iHmass, iP, *SatL, *SatV);
3434 CoolPropDbl dhV_dp = SatV->calc_first_saturation_deriv(iHmass, iP, *SatL, *SatV);
3435 CoolPropDbl dxdp_h = (Q() * dhV_dp + (1 - Q()) * dhL_dp) / (SatL->hmass() - SatV->hmass());
3436 CoolPropDbl dvdp_h = dvL_dp + dxdp_h * (1 / SatV->rhomass() - 1 / SatL->rhomass()) + Q() * (dvV_dp - dvL_dp);
3437 return -POW2(rhomass()) * dvdp_h;
3438 } else {
3439 throw ValueError("These inputs are not supported to calc_first_two_phase_deriv");
3440 }
3441}
3442CoolPropDbl HelmholtzEOSMixtureBackend::calc_first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, CoolPropDbl x_end) {
3443 // Note: If you need all three values (drho_dh__p, drho_dp__h and rho_spline),
3444 // you should calculate drho_dp__h first to avoid duplicate calculations.
3445
3446 bool drho_dh__p = false;
3447 bool drho_dp__h = false;
3448 bool rho_spline = false;
3449
3450 if (Of == iDmolar && Wrt == iHmolar && Constant == iP) {
3451 drho_dh__p = true;
3452 if (_drho_spline_dh__constp) return _drho_spline_dh__constp;
3453 } else if (Of == iDmass && Wrt == iHmass && Constant == iP) {
3454 return first_two_phase_deriv_splined(iDmolar, iHmolar, iP, x_end) * POW2(molar_mass());
3455 } else if (Of == iDmolar && Wrt == iP && Constant == iHmolar) {
3456 drho_dp__h = true;
3457 if (_drho_spline_dp__consth) return _drho_spline_dp__consth;
3458 } else if (Of == iDmass && Wrt == iP && Constant == iHmass) {
3459 return first_two_phase_deriv_splined(iDmolar, iP, iHmolar, x_end) * molar_mass();
3460 }
3461 // Add the special case for the splined density
3462 else if (Of == iDmolar && Wrt == iDmolar && Constant == iDmolar) {
3463 rho_spline = true;
3464 if (_rho_spline) return _rho_spline;
3465 } else if (Of == iDmass && Wrt == iDmass && Constant == iDmass) {
3466 return first_two_phase_deriv_splined(iDmolar, iDmolar, iDmolar, x_end) * molar_mass();
3467 } else {
3468 throw ValueError("These inputs are not supported to calc_first_two_phase_deriv");
3469 }
3470
3471 if (!this->SatL || !this->SatV) throw ValueError(format("The saturation properties are needed for calc_first_two_phase_deriv_splined"));
3472 if (_Q > x_end) {
3473 throw ValueError(format("Q [%g] is greater than x_end [%Lg]", _Q, x_end).c_str());
3474 }
3475 if (_phase != iphase_twophase) {
3476 throw ValueError(format("state is not two-phase"));
3477 }
3478
3479 shared_ptr<HelmholtzEOSMixtureBackend> Liq(new HelmholtzEOSMixtureBackend(this->get_components())),
3480 End(new HelmholtzEOSMixtureBackend(this->get_components()));
3481
3482 Liq->specify_phase(iphase_liquid);
3483 Liq->_Q = -1;
3484 Liq->update_DmolarT_direct(SatL->rhomolar(), SatL->T());
3485 End->update(QT_INPUTS, x_end, SatL->T());
3486
3487 CoolPropDbl Delta = Q() * (SatV->keyed_output(iHmolar) - SatL->keyed_output(iHmolar));
3488 CoolPropDbl Delta_end = End->keyed_output(iHmolar) - SatL->keyed_output(iHmolar);
3489
3490 // At the end of the zone to which spline is applied
3491 CoolPropDbl drho_dh_end = End->calc_first_two_phase_deriv(iDmolar, iHmolar, iP);
3492 CoolPropDbl rho_end = End->keyed_output(iDmolar);
3493
3494 // Faking single-phase
3495 CoolPropDbl rho_liq = Liq->keyed_output(iDmolar);
3496 CoolPropDbl drho_dh_liq__constp = Liq->first_partial_deriv(iDmolar, iHmolar, iP);
3497
3498 // Spline coordinates a, b, c, d
3499 CoolPropDbl Abracket = (2 * rho_liq - 2 * rho_end + Delta_end * (drho_dh_liq__constp + drho_dh_end));
3500 CoolPropDbl a = 1 / POW3(Delta_end) * Abracket;
3501 CoolPropDbl b = 3 / POW2(Delta_end) * (-rho_liq + rho_end) - 1 / Delta_end * (drho_dh_end + 2 * drho_dh_liq__constp);
3502 CoolPropDbl c = drho_dh_liq__constp;
3503 CoolPropDbl d = rho_liq;
3504
3505 // Either the spline value or drho/dh|p can be directly evaluated now
3506 _rho_spline = a * POW3(Delta) + b * POW2(Delta) + c * Delta + d;
3507 _drho_spline_dh__constp = 3 * a * POW2(Delta) + 2 * b * Delta + c;
3508 if (rho_spline) return _rho_spline;
3509 if (drho_dh__p) return _drho_spline_dh__constp;
3510
3511 // It's drho/dp|h
3512 // ... calculate some more things
3513
3514 // Derivatives *along* the saturation curve using the special internal method
3515 CoolPropDbl dhL_dp_sat = SatL->calc_first_saturation_deriv(iHmolar, iP, *SatL, *SatV);
3516 CoolPropDbl dhV_dp_sat = SatV->calc_first_saturation_deriv(iHmolar, iP, *SatL, *SatV);
3517 CoolPropDbl drhoL_dp_sat = SatL->calc_first_saturation_deriv(iDmolar, iP, *SatL, *SatV);
3518 CoolPropDbl drhoV_dp_sat = SatV->calc_first_saturation_deriv(iDmolar, iP, *SatL, *SatV);
3519 CoolPropDbl rhoV = SatV->keyed_output(iDmolar);
3520 CoolPropDbl rhoL = SatL->keyed_output(iDmolar);
3521 CoolPropDbl drho_dp_end = POW2(End->keyed_output(iDmolar)) * (x_end / POW2(rhoV) * drhoV_dp_sat + (1 - x_end) / POW2(rhoL) * drhoL_dp_sat);
3522
3523 // Faking single-phase
3524 //CoolPropDbl drho_dp__consth_liq = Liq->first_partial_deriv(iDmolar, iP, iHmolar);
3525 CoolPropDbl d2rhodhdp_liq = Liq->second_partial_deriv(iDmolar, iHmolar, iP, iP, iHmolar); // ?
3526
3527 // Derivatives at the end point
3528 // CoolPropDbl drho_dp__consth_end = End->calc_first_two_phase_deriv(iDmolar, iP, iHmolar);
3529 CoolPropDbl d2rhodhdp_end = End->calc_second_two_phase_deriv(iDmolar, iHmolar, iP, iP, iHmolar);
3530
3531 // Reminder:
3532 // Delta = Q()*(hV-hL) = h-hL
3533 // Delta_end = x_end*(hV-hL);
3534 CoolPropDbl d_Delta_dp__consth = -dhL_dp_sat;
3535 CoolPropDbl d_Delta_end_dp__consth = x_end * (dhV_dp_sat - dhL_dp_sat);
3536
3537 // First pressure derivative at constant h of the coefficients a,b,c,d
3538 // CoolPropDbl Abracket = (2*rho_liq - 2*rho_end + Delta_end * (drho_dh_liq__constp + drho_dh_end));
3539 CoolPropDbl d_Abracket_dp_consth = (2 * drhoL_dp_sat - 2 * drho_dp_end + Delta_end * (d2rhodhdp_liq + d2rhodhdp_end)
3540 + d_Delta_end_dp__consth * (drho_dh_liq__constp + drho_dh_end));
3541 CoolPropDbl da_dp = 1 / POW3(Delta_end) * d_Abracket_dp_consth + Abracket * (-3 / POW4(Delta_end) * d_Delta_end_dp__consth);
3542 CoolPropDbl db_dp = -6 / POW3(Delta_end) * d_Delta_end_dp__consth * (rho_end - rho_liq) + 3 / POW2(Delta_end) * (drho_dp_end - drhoL_dp_sat)
3543 + (1 / POW2(Delta_end) * d_Delta_end_dp__consth) * (drho_dh_end + 2 * drho_dh_liq__constp)
3544 - (1 / Delta_end) * (d2rhodhdp_end + 2 * d2rhodhdp_liq);
3545 CoolPropDbl dc_dp = d2rhodhdp_liq;
3546 CoolPropDbl dd_dp = drhoL_dp_sat;
3547
3548 _drho_spline_dp__consth =
3549 (3 * a * POW2(Delta) + 2 * b * Delta + c) * d_Delta_dp__consth + POW3(Delta) * da_dp + POW2(Delta) * db_dp + Delta * dc_dp + dd_dp;
3550 if (drho_dp__h) return _drho_spline_dp__consth;
3551
3552 throw ValueError("Something went wrong in HelmholtzEOSMixtureBackend::calc_first_two_phase_deriv_splined");
3553 return _HUGE;
3554}
3555
3556CoolProp::CriticalState HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0) {
3557 class Resid : public FuncWrapperND
3558 {
3559 public:
3560 HelmholtzEOSMixtureBackend& HEOS;
3561 double L1, M1;
3562 Eigen::MatrixXd Lstar, Mstar;
3563 Resid(HelmholtzEOSMixtureBackend& HEOS) : HEOS(HEOS), L1(_HUGE), M1(_HUGE){};
3564 std::vector<double> call(const std::vector<double>& tau_delta) {
3565 double rhomolar = tau_delta[1] * HEOS.rhomolar_reducing();
3566 double T = HEOS.T_reducing() / tau_delta[0];
3567 HEOS.update(DmolarT_INPUTS, rhomolar, T);
3568 Lstar = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT);
3569 Mstar = MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT, Lstar);
3570 std::vector<double> o(2);
3571 o[0] = Lstar.determinant();
3572 o[1] = Mstar.determinant();
3573 return o;
3574 };
3575 std::vector<std::vector<double>> Jacobian(const std::vector<double>& x) {
3576 std::size_t N = x.size();
3577 std::vector<std::vector<double>> J(N, std::vector<double>(N, 0));
3578 Eigen::MatrixXd adjL = adjugate(Lstar), adjM = adjugate(Mstar), dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau),
3579 dLdDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta),
3580 dMdTau = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iTau, Lstar, dLdTau),
3581 dMdDelta = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iDelta, Lstar, dLdDelta);
3582
3583 J[0][0] = (adjL * dLdTau).trace();
3584 J[0][1] = (adjL * dLdDelta).trace();
3585 J[1][0] = (adjM * dMdTau).trace();
3586 J[1][1] = (adjM * dMdDelta).trace();
3587 return J;
3588 }
3590 std::vector<std::vector<double>> numerical_Jacobian(const std::vector<double>& x) {
3591 std::size_t N = x.size();
3592 std::vector<double> rplus, rminus, xp;
3593 std::vector<std::vector<double>> J(N, std::vector<double>(N, 0));
3594
3595 double epsilon = 0.0001;
3596
3597 // Build the Jacobian by column
3598 for (std::size_t i = 0; i < N; ++i) {
3599 xp = x;
3600 xp[i] += epsilon;
3601 rplus = call(xp);
3602 xp[i] -= 2 * epsilon;
3603 rminus = call(xp);
3604
3605 for (std::size_t j = 0; j < N; ++j) {
3606 J[j][i] = (rplus[j] - rminus[j]) / (2 * epsilon);
3607 }
3608 }
3609 std::cout << J[0][0] << " " << J[0][1] << std::endl;
3610 std::cout << J[1][0] << " " << J[1][1] << std::endl;
3611 return J;
3612 };
3613 };
3614 Resid resid(*this);
3615 std::vector<double> x, tau_delta(2);
3616 tau_delta[0] = T_reducing() / T0;
3617 tau_delta[1] = rho0 / rhomolar_reducing();
3618 x = NDNewtonRaphson_Jacobian(&resid, tau_delta, 1e-10, 100);
3619 _critical.T = T_reducing() / x[0];
3620 _critical.rhomolar = x[1] * rhomolar_reducing();
3621 _critical.p = calc_pressure_nocache(_critical.T, _critical.rhomolar);
3622
3623 CriticalState critical;
3624 critical.T = _critical.T;
3625 critical.p = _critical.p;
3626 critical.rhomolar = _critical.rhomolar;
3627 if (_critical.p < 0) {
3628 critical.stable = false;
3629 } else {
3630 if (get_config_bool(ASSUME_CRITICAL_POINT_STABLE)) {
3631 critical.stable = true;
3632 } else {
3633 // Otherwise we try to check stability with TPD-based analysis
3634 StabilityRoutines::StabilityEvaluationClass stability_tester(*this);
3635 critical.stable = stability_tester.is_stable();
3636 }
3637 }
3638 return critical;
3639}
3640
3644class OneDimObjective : public FuncWrapper1DWithTwoDerivs
3645{
3646 public:
3647 CoolProp::HelmholtzEOSMixtureBackend& HEOS;
3648 const double delta;
3649 double _call, _deriv, _second_deriv;
3650 OneDimObjective(HelmholtzEOSMixtureBackend& HEOS, double delta0) : HEOS(HEOS), delta(delta0), _call(_HUGE), _deriv(_HUGE), _second_deriv(_HUGE){};
3651 double call(double tau) {
3652 double rhomolar = HEOS.rhomolar_reducing() * delta, T = HEOS.T_reducing() / tau;
3653 HEOS.update_DmolarT_direct(rhomolar, T);
3654 _call = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT).determinant();
3655 return _call;
3656 }
3657 double deriv(double tau) {
3658 Eigen::MatrixXd adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT)),
3659 dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau);
3660 _deriv = (adjL * dLdTau).trace();
3661 return _deriv;
3662 };
3663 double second_deriv(double tau) {
3664 Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT),
3665 dLstardTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau),
3666 d2LstardTau2 = MixtureDerivatives::d2Lstar_dX2(HEOS, XN_INDEPENDENT, iTau, iTau), adjL = adjugate(Lstar),
3667 dadjLstardTau = adjugate_derivative(Lstar, dLstardTau);
3668 _second_deriv = (dLstardTau * dadjLstardTau + adjL * d2LstardTau2).trace();
3669 return _second_deriv;
3670 };
3671};
3672
3675class L0CurveTracer : public FuncWrapper1DWithDeriv
3676{
3677 public:
3678 CoolProp::HelmholtzEOSMixtureBackend& HEOS;
3679 double delta, tau,
3680 M1_last,
3681 theta_last,
3682 R_tau,
3683 R_delta,
3684 R_tau_tracer,
3685 R_delta_tracer;
3686 std::vector<CoolProp::CriticalState> critical_points;
3687 int N_critical_points;
3688 Eigen::MatrixXd Lstar, adjLstar, dLstardTau, d2LstardTau2, dLstardDelta;
3689 SpinodalData spinodal_values;
3690 bool
3691 find_critical_points;
3692 L0CurveTracer(HelmholtzEOSMixtureBackend& HEOS, double tau0, double delta0)
3693 : HEOS(HEOS), delta(delta0), tau(tau0), M1_last(_HUGE), N_critical_points(0), find_critical_points(true) {
3694 R_delta_tracer = 0.1;
3695 R_delta = R_delta_tracer;
3696 R_tau_tracer = 0.1;
3697 R_tau = R_tau_tracer;
3698 };
3699 /***
3700 \brief Update values for tau, delta
3701 @param theta The angle
3702 @param tau The old value of tau
3703 @param delta The old value of delta
3704 @param tau_new The new value of tau
3705 @param delta_new The new value of delta
3706 */
3707 void get_tau_delta(const double theta, const double tau, const double delta, double& tau_new, double& delta_new) {
3708 tau_new = tau + R_tau * cos(theta);
3709 delta_new = delta + R_delta * sin(theta);
3710 };
3711 /***
3712 \brief Calculate the value of L1
3713 @param theta The angle
3714 */
3715 double call(double theta) {
3716 double tau_new, delta_new;
3717 this->get_tau_delta(theta, tau, delta, tau_new, delta_new);
3718 double rhomolar = HEOS.rhomolar_reducing() * delta_new, T = HEOS.T_reducing() / tau_new;
3719 HEOS.update_DmolarT_direct(rhomolar, T);
3720 Lstar = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT);
3721 adjLstar = adjugate(Lstar);
3722 dLstardTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau);
3723 dLstardDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta);
3724 double L1 = Lstar.determinant();
3725 return L1;
3726 };
3727 /***
3728 \brief Calculate the first partial derivative of L1 with respect to the angle
3729 @param theta The angle
3730 */
3731 double deriv(double theta) {
3732 double dL1_dtau = (adjLstar * dLstardTau).trace(), dL1_ddelta = (adjLstar * dLstardDelta).trace();
3733 return -R_tau * sin(theta) * dL1_dtau + R_delta * cos(theta) * dL1_ddelta;
3734 };
3735
3736 void trace() {
3737 bool debug = (get_debug_level() > 0) | false;
3738 double theta;
3739 for (int i = 0; i < 300; ++i) {
3740 if (i == 0) {
3741 // In the first iteration, search all angles in the positive delta direction using a
3742 // bounded solver with a very small radius in order to not hit other L1*=0 contours
3743 // that are in the vicinity
3744 R_tau = 0.001;
3745 R_delta = 0.001;
3746 theta = Brent(this, 0, M_PI, DBL_EPSILON, 1e-10, 100);
3747 } else {
3748 // In subsequent iterations, you already have an excellent guess for the direction to
3749 // be searching in, use Newton's method to refine the solution since we also
3750 // have an analytic solution for the derivative
3751 R_tau = R_tau_tracer;
3752 R_delta = R_delta_tracer;
3753 try {
3754 theta = Newton(this, theta_last, 1e-10, 100);
3755 } catch (...) {
3756 if (N_critical_points > 0 && delta > 1.5 * critical_points[0].rhomolar / HEOS.rhomolar_reducing()) {
3757 // Stopping the search - probably we have a kink in the L1*=0 contour
3758 // caused by poor low-temperature behavior
3759 continue;
3760 }
3761 }
3762
3763 // If the solver takes a U-turn, going in the opposite direction of travel
3764 // this is not a good thing, and force a Brent's method solver call to make sure we keep
3765 // tracing in the same direction
3766 if (std::abs(angle_difference(theta, theta_last)) > M_PI / 2.0) {
3767 // We have found at least one critical point and we are now well above the density
3768 // associated with the first critical point
3769 if (N_critical_points > 0 && delta > 1.2 * critical_points[0].rhomolar / HEOS.rhomolar_reducing()) {
3770 // Stopping the search - probably we have a kink in the L1*=0 contour
3771 // caused by poor low-temperature behavior
3772 continue;
3773 } else {
3774 theta = Brent(this, theta_last - M_PI / 3.5, theta_last + M_PI / 3.5, DBL_EPSILON, 1e-10, 100);
3775 }
3776 }
3777 }
3778
3779 // Calculate the second criticality condition
3780 double M1 = MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT, Lstar).determinant();
3781 double p_MPa = HEOS.p() / 1e6;
3782
3783 // Calculate the new tau and delta at the new point
3784 double tau_new, delta_new;
3785 this->get_tau_delta(theta, tau, delta, tau_new, delta_new);
3786
3787 // Stop if bounds are exceeded
3788 if (p_MPa > 500 || HEOS.get_critical_is_terminated(delta_new, tau_new)) {
3789 break;
3790 }
3791
3792 // If the sign of M1 and the previous value of M1 have different signs, it means that
3793 // you have bracketed a critical point, run the full critical point solver to
3794 // find the critical point and store it
3795 // Only enabled if find_critical_points is true (the default)
3796 if (i > 0 && M1 * M1_last < 0 && find_critical_points) {
3797 double rhomolar = HEOS.rhomolar_reducing() * (delta + delta_new) / 2.0, T = HEOS.T_reducing() / ((tau + tau_new) / 2.0);
3798 CoolProp::CriticalState crit = HEOS.calc_critical_point(rhomolar, T);
3799 critical_points.push_back(crit);
3800 N_critical_points++;
3801 if (debug) {
3802 std::cout << HEOS.get_mole_fractions()[0] << " " << crit.rhomolar << " " << crit.T << " " << p_MPa << std::endl;
3803 }
3804 }
3805
3806 // Update the storage values
3807 this->tau = tau_new;
3808 this->delta = delta_new;
3809 this->M1_last = M1;
3810 this->theta_last = theta;
3811
3812 this->spinodal_values.tau.push_back(tau_new);
3813 this->spinodal_values.delta.push_back(delta_new);
3814 this->spinodal_values.M1.push_back(M1);
3815 }
3816 };
3817};
3818
3819void HelmholtzEOSMixtureBackend::calc_criticality_contour_values(double& L1star, double& M1star) {
3820 Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(*this, XN_INDEPENDENT);
3821 Eigen::MatrixXd Mstar = MixtureDerivatives::Mstar(*this, XN_INDEPENDENT, Lstar);
3822 L1star = Lstar.determinant();
3823 M1star = Mstar.determinant();
3824};
3825
3826void HelmholtzEOSMixtureBackend::get_critical_point_search_radii(double& R_delta, double& R_tau) {
3827 R_delta = 0.025;
3828 R_tau = 0.1;
3829}
3830std::vector<CoolProp::CriticalState> HelmholtzEOSMixtureBackend::_calc_all_critical_points(bool find_critical_points) {
3831 // Populate the temporary class used to calculate the critical point(s)
3832 add_critical_state();
3833 if (get_debug_level() > 10) {
3834 rapidjson::Document doc;
3835 doc.SetObject();
3836 rapidjson::Value& val = doc;
3837 std::vector<std::vector<DepartureFunctionPointer>>& mat = critical_state->residual_helmholtz->Excess.DepartureFunctionMatrix;
3838 if (mat.size() > 0) {
3839 mat[0][1]->phi.to_json(val, doc);
3840 std::cout << cpjson::to_string(doc);
3841 }
3842 }
3843 critical_state->set_mole_fractions(this->get_mole_fractions_ref());
3844
3845 // Specify state to be something homogeneous to shortcut phase evaluation
3846 critical_state->specify_phase(iphase_gas);
3847
3848 double delta0 = _HUGE, tau0 = _HUGE;
3849 critical_state->get_critical_point_starting_values(delta0, tau0);
3850
3851 OneDimObjective resid_L0(*critical_state, delta0);
3852
3853 // If the derivative of L1star with respect to tau is positive,
3854 // tau needs to be increased such that we sit on the other
3855 // side of the d(L1star)/dtau = 0 contour
3856 resid_L0.call(tau0);
3857 int bump_count = 0;
3858 while (resid_L0.deriv(tau0) > 0 && bump_count < 3) {
3859 tau0 *= 1.1;
3860 bump_count++;
3861 }
3862 double tau_L0 = Halley(resid_L0, tau0, 1e-10, 100);
3863 //double T0 = T_reducing()/tau_L0;
3864 //double rho0 = delta0*rhomolar_reducing();
3865
3866 L0CurveTracer tracer(*critical_state, tau_L0, delta0);
3867 tracer.find_critical_points = find_critical_points;
3868
3869 double R_delta = 0, R_tau = 0;
3870 critical_state->get_critical_point_search_radii(R_delta, R_tau);
3871 tracer.R_delta_tracer = R_delta;
3872 tracer.R_tau_tracer = R_tau;
3873 tracer.trace();
3874
3875 this->spinodal_values = tracer.spinodal_values;
3876
3877 return tracer.critical_points;
3878}
3879
3880double HelmholtzEOSMixtureBackend::calc_tangent_plane_distance(const double T, const double p, const std::vector<double>& w,
3881 const double rhomolar_guess) {
3882
3883 const std::vector<CoolPropDbl>& z = this->get_mole_fractions_ref();
3884 if (w.size() != z.size()) {
3885 throw ValueError(format("Trial composition vector size [%d] is not the same as bulk composition [%d]", w.size(), z.size()));
3886 }
3887 add_TPD_state();
3888 TPD_state->set_mole_fractions(w);
3889
3890 CoolPropDbl rho = TPD_state->solver_rho_Tp_global(T, p, 0.9 / TPD_state->SRK_covolume());
3891 TPD_state->update_DmolarT_direct(rho, T);
3892
3893 double summer = 0;
3894 for (std::size_t i = 0; i < w.size(); ++i) {
3895 summer +=
3896 w[i] * (log(MixtureDerivatives::fugacity_i(*TPD_state, i, XN_DEPENDENT)) - log(MixtureDerivatives::fugacity_i(*this, i, XN_DEPENDENT)));
3897 }
3898 return summer;
3899}
3900
3901void HelmholtzEOSMixtureBackend::calc_build_spinodal() {
3902 // Ok, we are faking a little bit here, hijacking the code for critical points, but skipping evaluation of critical points
3903 bool find_critical_points = false;
3904 _calc_all_critical_points(find_critical_points);
3905}
3906
3907void HelmholtzEOSMixtureBackend::set_reference_stateS(const std::string& reference_state) {
3908 for (std::size_t i = 0; i < components.size(); ++i) {
3909 CoolProp::HelmholtzEOSMixtureBackend HEOS(std::vector<CoolPropFluid>(1, components[i]));
3910 if (!reference_state.compare("IIR")) {
3911 if (HEOS.Ttriple() > 273.15) {
3912 throw ValueError(format("Cannot use IIR reference state; Ttriple [%Lg] is greater than 273.15 K", HEOS.Ttriple()));
3913 }
3914 HEOS.update(QT_INPUTS, 0, 273.15);
3915
3916 // Get current values for the enthalpy and entropy
3917 double deltah = HEOS.hmass() - 200000; // offset from 200000 J/kg enthalpy
3918 double deltas = HEOS.smass() - 1000; // offset from 1000 J/kg/K entropy
3919 double delta_a1 = deltas / (HEOS.gas_constant() / HEOS.molar_mass());
3920 double delta_a2 = -deltah / (HEOS.gas_constant() / HEOS.molar_mass() * HEOS.get_reducing_state().T);
3921 // Change the value in the library for the given fluid
3922 set_fluid_enthalpy_entropy_offset(components[i], delta_a1, delta_a2, "IIR");
3923 if (get_debug_level() > 0) {
3924 std::cout << format("set offsets to %0.15g and %0.15g\n", delta_a1, delta_a2);
3925 }
3926 } else if (!reference_state.compare("ASHRAE")) {
3927 if (HEOS.Ttriple() > 233.15) {
3928 throw ValueError(format("Cannot use ASHRAE reference state; Ttriple [%Lg] is greater than than 233.15 K", HEOS.Ttriple()));
3929 }
3930 HEOS.update(QT_INPUTS, 0, 233.15);
3931
3932 // Get current values for the enthalpy and entropy
3933 double deltah = HEOS.hmass() - 0; // offset from 0 J/kg enthalpy
3934 double deltas = HEOS.smass() - 0; // offset from 0 J/kg/K entropy
3935 double delta_a1 = deltas / (HEOS.gas_constant() / HEOS.molar_mass());
3936 double delta_a2 = -deltah / (HEOS.gas_constant() / HEOS.molar_mass() * HEOS.get_reducing_state().T);
3937 // Change the value in the library for the given fluid
3938 set_fluid_enthalpy_entropy_offset(components[i], delta_a1, delta_a2, "ASHRAE");
3939 if (get_debug_level() > 0) {
3940 std::cout << format("set offsets to %0.15g and %0.15g\n", delta_a1, delta_a2);
3941 }
3942 } else if (!reference_state.compare("NBP")) {
3943 if (HEOS.p_triple() > 101325) {
3944 throw ValueError(format("Cannot use NBP reference state; p_triple [%Lg Pa] is greater than than 101325 Pa", HEOS.p_triple()));
3945 }
3946 // Saturated liquid boiling point at 1 atmosphere
3947 HEOS.update(PQ_INPUTS, 101325, 0);
3948
3949 double deltah = HEOS.hmass() - 0; // offset from 0 kJ/kg enthalpy
3950 double deltas = HEOS.smass() - 0; // offset from 0 kJ/kg/K entropy
3951 double delta_a1 = deltas / (HEOS.gas_constant() / HEOS.molar_mass());
3952 double delta_a2 = -deltah / (HEOS.gas_constant() / HEOS.molar_mass() * HEOS.get_reducing_state().T);
3953 // Change the value in the library for the given fluid
3954 set_fluid_enthalpy_entropy_offset(components[i], delta_a1, delta_a2, "NBP");
3955 if (get_debug_level() > 0) {
3956 std::cout << format("set offsets to %0.15g and %0.15g\n", delta_a1, delta_a2);
3957 }
3958 } else if (!reference_state.compare("DEF")) {
3959 set_fluid_enthalpy_entropy_offset(components[i], 0, 0, "DEF");
3960 } else if (!reference_state.compare("RESET")) {
3961 set_fluid_enthalpy_entropy_offset(components[i], 0, 0, "RESET");
3962 } else {
3963 throw ValueError(format("reference state string is invalid: [%s]", reference_state.c_str()));
3964 }
3965 }
3966}
3967
3973void HelmholtzEOSMixtureBackend::set_reference_stateD(double T, double rhomolar, double hmolar0, double smolar0) {
3974 for (std::size_t i = 0; i < components.size(); ++i) {
3975 CoolProp::HelmholtzEOSMixtureBackend HEOS(std::vector<CoolPropFluid>(1, components[i]));
3976
3977 HEOS.update(DmolarT_INPUTS, rhomolar, T);
3978
3979 // Get current values for the enthalpy and entropy
3980 double deltah = HEOS.hmolar() - hmolar0; // offset from specified enthalpy in J/mol
3981 double deltas = HEOS.smolar() - smolar0; // offset from specified entropy in J/mol/K
3982 double delta_a1 = deltas / (HEOS.gas_constant());
3983 double delta_a2 = -deltah / (HEOS.gas_constant() * HEOS.get_reducing_state().T);
3984 set_fluid_enthalpy_entropy_offset(components[i], delta_a1, delta_a2, "custom");
3985 }
3986}
3987
3988void HelmholtzEOSMixtureBackend::set_fluid_enthalpy_entropy_offset(CoolPropFluid& component, double delta_a1, double delta_a2,
3989 const std::string& ref) {
3990 component.EOS().alpha0.EnthalpyEntropyOffset.set(delta_a1, delta_a2, ref);
3991
3992 shared_ptr<CoolProp::HelmholtzEOSBackend> HEOS(new CoolProp::HelmholtzEOSBackend(component));
3993 HEOS->specify_phase(iphase_gas); // Something homogeneous;
3994 // Calculate the new enthalpy and entropy values
3995 HEOS->update(DmolarT_INPUTS, component.EOS().hs_anchor.rhomolar, component.EOS().hs_anchor.T);
3996 component.EOS().hs_anchor.hmolar = HEOS->hmolar();
3997 component.EOS().hs_anchor.smolar = HEOS->smolar();
3998
3999 double f = (HEOS->name() == "Water" || HEOS->name() == "CarbonDioxide") ? 1.00001 : 1.0;
4000
4001 // Calculate the new enthalpy and entropy values at the reducing state
4002 HEOS->update(DmolarT_INPUTS, component.EOS().reduce.rhomolar * f, component.EOS().reduce.T * f);
4003 component.EOS().reduce.hmolar = HEOS->hmolar();
4004 component.EOS().reduce.smolar = HEOS->smolar();
4005
4006 // Calculate the new enthalpy and entropy values at the critical state
4007 HEOS->update(DmolarT_INPUTS, component.crit.rhomolar * f, component.crit.T * f);
4008 component.crit.hmolar = HEOS->hmolar();
4009 component.crit.smolar = HEOS->smolar();
4010
4011 // Calculate the new enthalpy and entropy values
4012 HEOS->update(DmolarT_INPUTS, component.triple_liquid.rhomolar, component.triple_liquid.T);
4013 component.triple_liquid.hmolar = HEOS->hmolar();
4014 component.triple_liquid.smolar = HEOS->smolar();
4015
4016 // Calculate the new enthalpy and entropy values
4017 HEOS->update(DmolarT_INPUTS, component.triple_vapor.rhomolar, component.triple_vapor.T);
4018 component.triple_vapor.hmolar = HEOS->hmolar();
4019 component.triple_vapor.smolar = HEOS->smolar();
4020
4021 if (!HEOS->is_pure()) {
4022 // Calculate the new enthalpy and entropy values
4023 HEOS->update(DmolarT_INPUTS, component.EOS().max_sat_T.rhomolar, component.EOS().max_sat_T.T);
4024 component.EOS().max_sat_T.hmolar = HEOS->hmolar();
4025 component.EOS().max_sat_T.smolar = HEOS->smolar();
4026 // Calculate the new enthalpy and entropy values
4027 HEOS->update(DmolarT_INPUTS, component.EOS().max_sat_p.rhomolar, component.EOS().max_sat_p.T);
4028 component.EOS().max_sat_p.hmolar = HEOS->hmolar();
4029 component.EOS().max_sat_p.smolar = HEOS->smolar();
4030 }
4031}
4032
4033} /* namespace CoolProp */
Generated on Mon Jun 2 2025 05:21:52 for CoolProp by doxygen 1.9.4