CoolProp 8.0.0
An open-source fluid property and humid air property database
HelmholtzEOSMixtureBackend.h
Go to the documentation of this file.
1
2#ifndef HELMHOLTZEOSMIXTUREBACKEND_H_
3#define HELMHOLTZEOSMIXTUREBACKEND_H_
4
5#include "CoolProp/AbstractState.h"
6#include "CoolProp/CoolPropFluid.h"
7#include "ReducingFunctions.h"
8#include "ExcessHEFunction.h"
9#include "CoolProp/numerics/Solvers.h"
10#include "CoolProp/fluids/PhaseEnvelope.h"
11#include "CoolProp/DataStructures.h"
12#include "CoolProp/Configuration.h"
13
14#include <vector>
15
16namespace CoolProp {
17
18class FlashRoutines;
19
20class ResidualHelmholtz;
21
22// This class contains the mole fractions for a given mixture.
23class MoleFractions
24{
25 private:
26 std::vector<CoolPropDbl> mole_fractions;
27 template <typename T>
28 bool verify_mole_fractions_set(T i) const {
29 if (i >= mole_fractions.size()) {
30 throw CoolProp::ValueError("mole fractions are not set for all components");
31 }
32 return true;
33 }
34
35 public:
36 template <typename T>
37 void resize(T N) {
38 return mole_fractions.resize(N);
39 }
40 std::size_t size() const {
41 return mole_fractions.size();
42 }
43 void clear() {
44 mole_fractions.clear();
45 }
46 // operator overloads
47 template <typename T>
48 MoleFractions& operator=(const std::vector<T>& values) {
49 mole_fractions = values;
50 return *this;
51 }
52 template <typename T>
53 CoolPropDbl operator[](T i) const {
54 verify_mole_fractions_set(i);
55 return mole_fractions[i];
56 }
57 operator std::vector<CoolPropDbl>&() {
58 return mole_fractions;
59 }
60};
61
62class HelmholtzEOSMixtureBackend : public AbstractState
63{
64
65 protected:
66 void pre_update(CoolProp::input_pairs& input_pair, CoolPropDbl& value1, CoolPropDbl& value2);
67 void post_update(bool optional_checks = true);
68 std::vector<shared_ptr<HelmholtzEOSMixtureBackend>>
69 linked_states;
70 shared_ptr<HelmholtzEOSMixtureBackend> transient_pure_state;
71 shared_ptr<HelmholtzEOSMixtureBackend> TPD_state;
72 shared_ptr<HelmholtzEOSMixtureBackend> critical_state;
78 shared_ptr<HelmholtzEOSMixtureBackend> viscosity_ecs_reference_state;
79 shared_ptr<HelmholtzEOSMixtureBackend> conductivity_ecs_reference_state;
81 virtual void add_TPD_state() {
82 if (TPD_state.get() == nullptr) {
83 bool sat_states = false;
84 TPD_state.reset(get_copy(sat_states));
85 linked_states.push_back(TPD_state);
86 }
87 };
89 virtual void add_critical_state() {
90 if (critical_state.get() == nullptr) {
91 bool sat_states = true;
92 critical_state.reset(get_copy(sat_states));
93 linked_states.push_back(critical_state);
94 }
95 };
97 virtual void add_transient_pure_state() {
98 if (transient_pure_state.get() == nullptr) {
99 bool sat_states = true;
100 transient_pure_state.reset(get_copy(sat_states));
101 linked_states.push_back(transient_pure_state);
102 }
103 };
104
105 std::vector<CoolPropFluid> components;
106 bool is_pure_or_pseudopure;
107 MoleFractions mole_fractions;
108 std::vector<CoolPropDbl> K,
109 lnK;
110
111 SimpleState _crit;
112 std::size_t N;
113
115 std::vector<CoolProp::CriticalState> _calc_all_critical_points(bool find_critical_points = true);
116
117 static void set_fluid_enthalpy_entropy_offset(CoolPropFluid& component, double delta_a1, double delta_a2, const std::string& ref);
118
119 public:
123 std::shared_ptr<EquationOfState::SuperAncillary_t> get_superanc();
124
133 void ensure_caloric_superancillaries();
134
135 HelmholtzEOSMixtureBackend();
136 HelmholtzEOSMixtureBackend(const std::vector<CoolPropFluid>& components, bool generate_SatL_and_SatV = true);
137 HelmholtzEOSMixtureBackend(const std::vector<std::string>& component_names, bool generate_SatL_and_SatV = true);
138 virtual HelmholtzEOSMixtureBackend* get_copy(bool generate_SatL_and_SatV = true);
139
140 // Copy over the reducing and departure terms to all linked states (recursively)
141 void sync_linked_states(const HelmholtzEOSMixtureBackend* const);
142
143 virtual ~HelmholtzEOSMixtureBackend() = default;
144 std::string backend_name() override {
145 return get_backend_string(HEOS_BACKEND_MIX);
146 }
147 shared_ptr<ReducingFunction> Reducing;
148 shared_ptr<ResidualHelmholtz> residual_helmholtz;
149 PhaseEnvelopeData PhaseEnvelope;
150 SimpleState hsat_max;
151 SsatSimpleState ssat_max;
152 SpinodalData spinodal_values;
153
154 bool clear() override {
155 // Clear the locally cached values for the derivatives of the Helmholtz energy
156 // in each component
157 for (auto& component : components) {
158 component.EOS().alphar.clear();
159 component.EOS().alpha0.clear();
160 }
161 return AbstractState::clear();
162 };
163
164 friend class
165 FlashRoutines; // Allows the static methods in the FlashRoutines class to have access to all the protected members and methods of this class
166 friend class
167 TransportRoutines; // Allows the static methods in the TransportRoutines class to have access to all the protected members and methods of this class
168 friend class
169 MixtureDerivatives; // Allows the static methods in the MixtureDerivatives class to have access to all the protected members and methods of this class
170 friend class
171 PhaseEnvelopeRoutines; // Allows the static methods in the PhaseEnvelopeRoutines class to have access to all the protected members and methods of this class
172 friend class
173 MixtureParameters; // Allows the static methods in the MixtureParameters class to have access to all the protected members and methods of this class
174 friend class
175 CorrespondingStatesTerm; // // Allows the methods in the CorrespondingStatesTerm class to have access to all the protected members and methods of this class
176
177 // Helmholtz EOS backend uses mole fractions
178 bool using_mole_fractions() override {
179 return true;
180 }
181 bool using_mass_fractions() override {
182 return false;
183 }
184 bool using_volu_fractions() override {
185 return false;
186 }
187 bool is_pure() {
188 return components.size() == 1 && !components[0].EOS().pseudo_pure;
189 }
190 bool has_melting_line() override {
191 return is_pure_or_pseudopure && components[0].ancillaries.melting_line.enabled();
192 };
193 CoolPropDbl calc_melting_line(int param, int given, CoolPropDbl value) override;
194
196 std::vector<std::pair<double, double>> get_melting_line_part_pranges() {
197 if (!is_pure_or_pseudopure) {
198 throw NotImplementedError("melting line not available for mixtures");
199 }
200 std::vector<std::pair<double, double>> out;
201 for (const auto& pr : components[0].ancillaries.melting_line.get_parts_pranges()) {
202 out.emplace_back(static_cast<double>(pr.first), static_cast<double>(pr.second));
203 }
204 return out;
205 }
206
208 std::string fluid_param_string(const std::string&) override;
209
211 void set_reference_stateS(const std::string& reference_state) override;
212
214 void set_reference_stateD(double T, double rhomolar, double hmolar0, double smolar0) override;
215
217 void set_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter, const double value) override;
219 double get_binary_interaction_double(const std::size_t i, const std::size_t j, const std::string& parameter) override;
221 //virtual std::string get_binary_interaction_string(const std::size_t &i, const std::size_t &j, const std::string &parameter);
223 void set_binary_interaction_string(const std::size_t i, const std::size_t j, const std::string& parameter, const std::string& value) override;
225 void apply_simple_mixing_rule(std::size_t i, std::size_t j, const std::string& model) override;
226
227 // Set the cubic alpha function's constants:
228 void set_cubic_alpha_C(const size_t i, const std::string& parameter, const double c1, const double c2, const double c3) override {
229 throw ValueError("set_cubic_alpha_C only defined for cubic backends");
230 };
231
232 // Set fluid parameter (currently the volume translation parameter for cubic)
233 void set_fluid_parameter_double(const size_t i, const std::string& parameter, const double value) override {
234 throw ValueError("set_fluid_parameter_double only defined for cubic backends");
235 };
236 double get_fluid_parameter_double(const size_t i, const std::string& parameter) override;
237
238 phases calc_phase() override {
239 return _phase;
240 };
241
246 void calc_specify_phase(phases phase_index) override {
247 imposed_phase_index = phase_index;
248 _phase = phase_index;
249 }
252 void calc_unspecify_phase() override {
253 imposed_phase_index = iphase_not_imposed;
254 }
255 CoolPropDbl calc_saturation_ancillary(parameters param, int Q, parameters given, double value) override;
256 void calc_ssat_max();
257 void calc_hsat_max();
258 CoolPropDbl calc_GWP20() override;
259 CoolPropDbl calc_GWP500() override;
260 CoolPropDbl calc_GWP100() override;
261 CoolPropDbl calc_ODP() override;
262
263 CoolPropDbl calc_first_saturation_deriv(parameters Of1, parameters Wrt1) override;
264 CoolPropDbl calc_first_saturation_deriv(parameters Of1, parameters Wrt1, HelmholtzEOSMixtureBackend& SatL, HelmholtzEOSMixtureBackend& SatV);
265 CoolPropDbl calc_second_saturation_deriv(parameters Of1, parameters Wrt1, parameters Wrt2) override;
266 CoolPropDbl calc_first_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant) override;
267 CoolPropDbl calc_second_two_phase_deriv(parameters Of, parameters Wrt1, parameters Constant1, parameters Wrt2, parameters Constant2) override;
268 CoolPropDbl calc_first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, CoolPropDbl x_end) override;
269
270 CriticalState calc_critical_point(double rho0, double T0);
272 std::vector<CoolProp::CriticalState> calc_all_critical_points() override {
273 bool find_critical_points = true;
274 return _calc_all_critical_points(find_critical_points);
275 };
276
277 virtual void get_critical_point_starting_values(double& delta0, double& tau0) {
278 delta0 = get_config_double(SPINODAL_MINIMUM_DELTA); // The value of delta where we start searching for crossing with Lstar=0 contour
279 tau0 = 0.66; // The value of tau where we start searching at delta=delta0
280 }
282 virtual void get_critical_point_search_radii(double& R_delta, double& R_tau);
284 virtual bool get_critical_is_terminated(double& delta, double& tau) {
285 return delta > 5 || tau > 5;
286 }
287
289 void calc_build_spinodal() override;
290
292 SpinodalData calc_get_spinodal_data() override {
293 return spinodal_values;
294 };
295
297 void calc_criticality_contour_values(double& L1star, double& M1star) override;
298
300 double calc_tangent_plane_distance(const double T, const double p, const std::vector<double>& w, const double rhomolar_guess) override;
301
303 void recalculate_singlephase_phase();
304
306 void calc_change_EOS(const std::size_t i, const std::string& EOS_name) override;
307
308 const CoolProp::SimpleState& calc_state(const std::string& state) override;
309
310 const double get_fluid_constant(std::size_t i, parameters param) const override {
311 const CoolPropFluid& fld = components[i];
312 switch (param) {
313 case iP_critical:
314 return fld.crit.p;
315 case iT_critical:
316 return fld.crit.T;
317 case iT_reducing:
318 return fld.EOS().reduce.T;
319 case irhomolar_reducing:
320 return fld.EOS().reduce.rhomolar;
321 case irhomolar_critical:
322 return fld.crit.rhomolar;
323 case iacentric_factor:
324 return fld.EOS().acentric;
325 case imolar_mass:
326 return fld.EOS().molar_mass;
327 case iT_triple:
328 return fld.EOS().sat_min_liquid.T;
329 case iP_triple:
330 return fld.EOS().sat_min_liquid.p;
331 case igas_constant:
332 return fld.EOS().R_u;
333 default:
334 throw ValueError(format("I don't know what to do with this fluid constant: %s", get_parameter_information(param, "short").c_str()));
335 }
336 }
337
338 const std::vector<CoolPropFluid>& get_components() const {
339 return components;
340 }
341 std::vector<CoolPropFluid>& get_components() {
342 return components;
343 }
344 std::vector<CoolPropDbl>& get_K() {
345 return K;
346 };
347 std::vector<CoolPropDbl>& get_lnK() {
348 return lnK;
349 };
350 HelmholtzEOSMixtureBackend& get_SatL() {
351 return *SatL;
352 };
353 HelmholtzEOSMixtureBackend& get_SatV() {
354 return *SatV;
355 };
356
357 std::vector<CoolPropDbl> calc_mole_fractions_liquid() override {
358 // SatL/SatV retain composition vectors from the most recent VLE flash
359 // (or phase-envelope build) — those are not meaningful when the
360 // current state is single-phase, and returning them silently surfaced
361 // as bug #2308. Require the current state to be two-phase before
362 // returning anything.
363 if (_phase != iphase_twophase) {
364 throw ValueError("mole_fractions_liquid is only defined in the two-phase region (current state is single-phase)");
365 }
366 return SatL->get_mole_fractions();
367 };
368 std::vector<CoolPropDbl> calc_mole_fractions_vapor() override {
369 if (_phase != iphase_twophase) {
370 throw ValueError("mole_fractions_vapor is only defined in the two-phase region (current state is single-phase)");
371 }
372 return SatV->get_mole_fractions();
373 };
374
375 const std::vector<CoolPropDbl> calc_mass_fractions() override;
376
377 const CoolProp::PhaseEnvelopeData& calc_phase_envelope_data() override {
378 return PhaseEnvelope;
379 };
380
382 void calc_conformal_state(const std::string& reference_fluid, CoolPropDbl& T, CoolPropDbl& rhomolar) override;
383
384 void resize(std::size_t N);
385 shared_ptr<HelmholtzEOSMixtureBackend> SatL, SatV;
386
392 void update(CoolProp::input_pairs input_pair, double value1, double value2) override;
393
397 void update_with_guesses(CoolProp::input_pairs input_pair, double Value1, double Value2, const GuessesStructure& guesses) override;
398
401 void update_internal(HelmholtzEOSMixtureBackend& HEOS);
402
408 void update_TP_guessrho(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rho_guess);
409 void update_DmolarT_direct(CoolPropDbl rhomolar, CoolPropDbl T);
410 void update_TDmolarP_unchecked(CoolPropDbl T, CoolPropDbl rhomolarL, CoolPropDbl p);
411 void update_QT_pure_superanc(CoolPropDbl Q, CoolPropDbl T) override;
412 void update_HmolarQ_with_guessT(CoolPropDbl hmolar, CoolPropDbl Q, CoolPropDbl Tguess);
413
419 virtual void set_components(const std::vector<CoolPropFluid>& components, bool generate_SatL_and_SatV = true);
420
423 void set_mixture_parameters();
424
429 void set_mole_fractions(const std::vector<CoolPropDbl>& mf) override;
430
431 const std::vector<CoolPropDbl>& get_mole_fractions() override {
432 return mole_fractions;
433 };
434 std::vector<CoolPropDbl>& get_mole_fractions_ref() {
435 return mole_fractions;
436 };
437 std::vector<double>& get_mole_fractions_doubleref() {
438 return mole_fractions;
439 }
440
445 void set_mass_fractions(const std::vector<CoolPropDbl>& mass_fractions) override;
446
447 void calc_ideal_curve(const std::string& type, std::vector<double>& T, std::vector<double>& p) override;
448
449 CoolPropDbl calc_molar_mass() override;
450 PhaseMolarMasses calc_phase_molar_masses() override;
451 CoolPropDbl calc_gas_constant() override;
452 CoolPropDbl calc_acentric_factor() override;
453
454 CoolPropDbl calc_Bvirial() override;
455 CoolPropDbl calc_Cvirial() override;
456 CoolPropDbl calc_dBvirial_dT() override;
457 CoolPropDbl calc_dCvirial_dT() override;
458
459 CoolPropDbl calc_pressure() override;
460 CoolPropDbl calc_cvmolar() override;
461 CoolPropDbl calc_cpmolar() override;
462 CoolPropDbl calc_gibbsmolar() override;
463 CoolPropDbl calc_gibbsmolar_residual() override {
464 return gas_constant() * _T * (alphar() + delta() * dalphar_dDelta());
465 }
466 CoolPropDbl calc_gibbsmolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar);
467
468 CoolPropDbl calc_helmholtzmolar() override;
469 CoolPropDbl calc_cpmolar_idealgas() override;
470 CoolPropDbl calc_pressure_nocache(CoolPropDbl T, CoolPropDbl rhomolar);
471 CoolPropDbl calc_smolar() override;
472 CoolPropDbl calc_smolar_residual() override {
473 return gas_constant() * (tau() * dalphar_dTau() - alphar());
474 }
475 CoolPropDbl calc_smolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar);
476
477 CoolPropDbl calc_hmolar() override;
478 CoolPropDbl calc_hmolar_residual() override {
479 return gas_constant() * _T * (tau() * dalphar_dTau() + delta() * dalphar_dDelta());
480 }
481 CoolPropDbl calc_hmolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar);
482
483 CoolPropDbl calc_umolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar);
484 CoolPropDbl calc_umolar() override;
485 CoolPropDbl calc_speed_sound() override;
486
487 void calc_excess_properties() override;
488
494 CoolPropDbl calc_phase_identification_parameter();
495 CoolPropDbl calc_fugacity(std::size_t i) override;
496 CoolPropDbl calc_fugacity_coefficient(std::size_t i) override;
497 CoolPropDbl calc_chemical_potential(std::size_t i) override;
498
500 CoolPropDbl calc_flame_hazard() override {
501 return components[0].environment.FH;
502 };
504 CoolPropDbl calc_health_hazard() override {
505 return components[0].environment.HH;
506 };
508 CoolPropDbl calc_physical_hazard() override {
509 return components[0].environment.PH;
510 };
511
513 CoolPropDbl calc_alphar() override;
515 CoolPropDbl calc_dalphar_dDelta() override;
517 CoolPropDbl calc_dalphar_dTau() override;
519 CoolPropDbl calc_d2alphar_dDelta2() override;
521 CoolPropDbl calc_d2alphar_dDelta_dTau() override;
523 CoolPropDbl calc_d2alphar_dTau2() override;
525 CoolPropDbl calc_d3alphar_dDelta3() override;
527 CoolPropDbl calc_d3alphar_dDelta2_dTau() override;
529 CoolPropDbl calc_d3alphar_dDelta_dTau2() override;
531 CoolPropDbl calc_d3alphar_dTau3() override;
533 CoolPropDbl calc_d4alphar_dDelta4() override;
535 CoolPropDbl calc_d4alphar_dDelta3_dTau() override;
537 CoolPropDbl calc_d4alphar_dDelta2_dTau2() override;
539 CoolPropDbl calc_d4alphar_dDelta_dTau3() override;
541 CoolPropDbl calc_d4alphar_dTau4() override;
542
543 CoolPropDbl calc_alpha0() override;
544 CoolPropDbl calc_dalpha0_dDelta() override;
545 CoolPropDbl calc_dalpha0_dTau() override;
546 CoolPropDbl calc_d2alpha0_dDelta2() override;
547 CoolPropDbl calc_d2alpha0_dDelta_dTau() override;
548 CoolPropDbl calc_d2alpha0_dTau2() override;
549 CoolPropDbl calc_d3alpha0_dDelta3() override;
550 CoolPropDbl calc_d3alpha0_dDelta2_dTau() override;
551 CoolPropDbl calc_d3alpha0_dDelta_dTau2() override;
552 CoolPropDbl calc_d3alpha0_dTau3() override;
553
554 CoolPropDbl calc_surface_tension() override;
555 CoolPropDbl calc_viscosity() override;
556 CoolPropDbl calc_viscosity_dilute();
557 CoolPropDbl calc_viscosity_background();
558 CoolPropDbl calc_viscosity_background(CoolPropDbl eta_dilute, CoolPropDbl& initial_density, CoolPropDbl& residual);
559 CoolPropDbl calc_conductivity() override;
560 CoolPropDbl calc_conductivity_background();
561
567 void calc_viscosity_contributions(CoolPropDbl& dilute, CoolPropDbl& initial_density, CoolPropDbl& residual, CoolPropDbl& critical) override;
573 void calc_conductivity_contributions(CoolPropDbl& dilute, CoolPropDbl& initial_density, CoolPropDbl& residual, CoolPropDbl& critical) override;
574
575 CoolPropDbl calc_saturated_liquid_keyed_output(parameters key) override;
576 CoolPropDbl calc_saturated_vapor_keyed_output(parameters key) override;
577
578 CoolPropDbl calc_Tmin() override;
579 CoolPropDbl calc_Tmax() override;
580 CoolPropDbl calc_pmax() override;
581 CoolPropDbl calc_Ttriple() override;
582 CoolPropDbl calc_p_triple() override;
583 CoolPropDbl calc_pmax_sat();
584 CoolPropDbl calc_Tmax_sat();
585 void calc_Tmin_sat(CoolPropDbl& Tmin_satL, CoolPropDbl& Tmin_satV);
586 void calc_pmin_sat(CoolPropDbl& pmin_satL, CoolPropDbl& pmin_satV);
587
588 CoolPropDbl calc_T_critical() override;
589 CoolPropDbl calc_p_critical() override;
590 CoolPropDbl calc_rhomolar_critical() override;
591
592 CoolPropDbl calc_T_reducing() override {
593 return get_reducing_state().T;
594 };
595 CoolPropDbl calc_rhomolar_reducing() override {
596 return get_reducing_state().rhomolar;
597 };
598 CoolPropDbl calc_p_reducing() override {
599 return get_reducing_state().p;
600 };
601
602 // Calculate the phase identification parameter of Venkatarathnam et al, Fluid Phase Equilibria
603 CoolPropDbl calc_PIP() override {
604 return 2
605 - rhomolar()
606 * (second_partial_deriv(iP, iDmolar, iT, iT, iDmolar) / first_partial_deriv(iP, iT, iDmolar)
607 - second_partial_deriv(iP, iDmolar, iT, iDmolar, iT) / first_partial_deriv(iP, iDmolar, iT));
608 };
609
610 std::string calc_name() override;
611 std::vector<std::string> calc_fluid_names() override;
612
613 void calc_all_alphar_deriv_cache(const std::vector<CoolPropDbl>& mole_fractions, const CoolPropDbl& tau, const CoolPropDbl& delta);
614 virtual CoolPropDbl calc_alphar_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl>& mole_fractions,
615 const CoolPropDbl& tau, const CoolPropDbl& delta);
616
641 CoolPropDbl calc_alpha0_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl>& mole_fractions, const CoolPropDbl& tau,
642 const CoolPropDbl& delta, const CoolPropDbl& Tr, const CoolPropDbl& rhor);
643
644 HelmholtzDerivatives calc_all_alpha0_derivs_nocache(const std::vector<CoolPropDbl>& mole_fractions, const CoolPropDbl& tau,
645 const CoolPropDbl& delta, const CoolPropDbl& Tr, const CoolPropDbl& rhor);
646
647 void calc_reducing_state() override;
648 virtual SimpleState calc_reducing_state_nocache(const std::vector<CoolPropDbl>& mole_fractions);
649
650 const CoolProp::SimpleState& get_reducing_state() override {
651 calc_reducing_state();
652 return _reducing;
653 };
654
655 void update_states() override;
656
657 CoolPropDbl calc_compressibility_factor() override {
658 return 1 + delta() * dalphar_dDelta();
659 };
660
661 void calc_phase_envelope(const std::string& type) override;
662
663 CoolPropDbl SRK_covolume();
664 virtual CoolPropDbl calc_rhomolar_max_bound();
665
666 // ***************************************************************
667 // ***************************************************************
668 // ************* PHASE DETERMINATION ROUTINES ******************
669 // ***************************************************************
670 // ***************************************************************
671 void T_phase_determination_pure_or_pseudopure(int other, CoolPropDbl value);
672 void p_phase_determination_pure_or_pseudopure(int other, CoolPropDbl value, bool& saturation_called);
673 void DmolarP_phase_determination();
674
675 // ***************************************************************
676 // ***************************************************************
677 // ******************* SOLVER ROUTINES *************************
678 // ***************************************************************
679 // ***************************************************************
680
681 virtual CoolPropDbl solver_rho_Tp(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rho_guess = -1);
682 virtual CoolPropDbl solver_rho_Tp_SRK(CoolPropDbl T, CoolPropDbl p, phases phase);
683 enum StationaryPointReturnFlag
684 {
685 ZERO_STATIONARY_POINTS,
686 ONE_STATIONARY_POINT_FOUND,
687 TWO_STATIONARY_POINTS_FOUND
688 };
689 virtual StationaryPointReturnFlag solver_dpdrho0_Tp(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomax, CoolPropDbl& light, CoolPropDbl& heavy);
690 virtual CoolPropDbl solver_rho_Tp_global(CoolPropDbl T, CoolPropDbl p, CoolPropDbl rhomax);
691};
692
693class CorrespondingStatesTerm
694{
695 public:
697 virtual HelmholtzDerivatives all(HelmholtzEOSMixtureBackend& HEOS, double tau, double delta, const std::vector<CoolPropDbl>& x,
698 bool cache_values = false) {
699 HelmholtzDerivatives summer;
700 std::size_t N = x.size();
701 for (std::size_t i = 0; i < N; ++i) {
702 HelmholtzDerivatives derivs = HEOS.components[i].EOS().alphar.all(tau, delta, cache_values);
703 summer = summer + derivs * x[i];
704 }
705 return summer;
706 }
707 CoolPropDbl dalphar_dxi(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
708 if (xN_flag == XN_INDEPENDENT) {
709 return HEOS.components[i].EOS().baser(HEOS.tau(), HEOS.delta());
710 } else if (xN_flag == XN_DEPENDENT) {
711 std::size_t N = x.size();
712 if (i == N - 1) return 0;
713 return HEOS.components[i].EOS().baser(HEOS.tau(), HEOS.delta()) - HEOS.components[N - 1].EOS().baser(HEOS.tau(), HEOS.delta());
714 } else {
715 throw ValueError(format("xN_flag is invalid"));
716 }
717 }
718 CoolPropDbl d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
719 if (xN_flag == XN_INDEPENDENT) {
720 return HEOS.components[i].EOS().dalphar_dTau(HEOS._tau, HEOS._delta);
721 } else if (xN_flag == XN_DEPENDENT) {
722 std::size_t N = x.size();
723 if (i == N - 1) return 0;
724 return HEOS.components[i].EOS().dalphar_dTau(HEOS._tau, HEOS._delta) - HEOS.components[N - 1].EOS().dalphar_dTau(HEOS._tau, HEOS._delta);
725 } else {
726 throw ValueError(format("xN_flag is invalid"));
727 }
728 }
729 CoolPropDbl d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
730 if (xN_flag == XN_INDEPENDENT) {
731 return HEOS.components[i].EOS().dalphar_dDelta(HEOS.tau(), HEOS.delta());
732 } else if (xN_flag == XN_DEPENDENT) {
733 std::size_t N = x.size();
734 if (i == N - 1) return 0;
735 return HEOS.components[i].EOS().dalphar_dDelta(HEOS.tau(), HEOS.delta())
736 - HEOS.components[N - 1].EOS().dalphar_dDelta(HEOS._tau, HEOS._delta);
737 } else {
738 throw ValueError(format("xN_flag is invalid"));
739 }
740 }
741 CoolPropDbl d3alphar_dxi_dDelta2(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
742 if (xN_flag == XN_INDEPENDENT) {
743 return HEOS.components[i].EOS().d2alphar_dDelta2(HEOS.tau(), HEOS.delta());
744 } else if (xN_flag == XN_DEPENDENT) {
745 std::size_t N = x.size();
746 if (i == N - 1) return 0;
747 return HEOS.components[i].EOS().d2alphar_dDelta2(HEOS.tau(), HEOS.delta())
748 - HEOS.components[N - 1].EOS().d2alphar_dDelta2(HEOS.tau(), HEOS.delta());
749 } else {
750 throw ValueError(format("xN_flag is invalid"));
751 }
752 }
753 CoolPropDbl d3alphar_dxi_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
754 if (xN_flag == XN_INDEPENDENT) {
755 return HEOS.components[i].EOS().d2alphar_dTau2(HEOS.tau(), HEOS.delta());
756 } else if (xN_flag == XN_DEPENDENT) {
757 std::size_t N = x.size();
758 if (i == N - 1) return 0;
759 return HEOS.components[i].EOS().d2alphar_dTau2(HEOS.tau(), HEOS.delta())
760 - HEOS.components[N - 1].EOS().d2alphar_dTau2(HEOS.tau(), HEOS.delta());
761 } else {
762 throw ValueError(format("xN_flag is invalid"));
763 }
764 }
765 CoolPropDbl d3alphar_dxi_dDelta_dTau(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
766 if (xN_flag == XN_INDEPENDENT) {
767 return HEOS.components[i].EOS().d2alphar_dDelta_dTau(HEOS.tau(), HEOS.delta());
768 } else if (xN_flag == XN_DEPENDENT) {
769 std::size_t N = x.size();
770 if (i == N - 1) return 0;
771 return HEOS.components[i].EOS().d2alphar_dDelta_dTau(HEOS.tau(), HEOS.delta())
772 - HEOS.components[N - 1].EOS().d2alphar_dDelta_dTau(HEOS.tau(), HEOS.delta());
773 } else {
774 throw ValueError(format("xN_flag is invalid"));
775 }
776 }
777
778 CoolPropDbl d2alphardxidxj(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
779 x_N_dependency_flag xN_flag) {
780 if (xN_flag == XN_INDEPENDENT) {
781 return 0;
782 } else if (xN_flag == XN_DEPENDENT) {
783 return 0;
784 } else {
785 throw ValueError(format("xN_flag is invalid"));
786 }
787 }
788
789 CoolPropDbl d4alphar_dxi_dDelta3(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
790 if (xN_flag == XN_INDEPENDENT) {
791 return HEOS.components[i].EOS().d3alphar_dDelta3(HEOS.tau(), HEOS.delta());
792 } else {
793 throw ValueError(format("xN_flag is invalid"));
794 }
795 }
796 CoolPropDbl d4alphar_dxi_dTau3(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
797 if (xN_flag == XN_INDEPENDENT) {
798 return HEOS.components[i].EOS().d3alphar_dTau3(HEOS.tau(), HEOS.delta());
799 } else {
800 throw ValueError(format("xN_flag is invalid"));
801 }
802 }
803 CoolPropDbl d4alphar_dxi_dDelta_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
804 if (xN_flag == XN_INDEPENDENT) {
805 return HEOS.components[i].EOS().d3alphar_dDelta_dTau2(HEOS.tau(), HEOS.delta());
806 } else {
807 throw ValueError(format("xN_flag is invalid"));
808 }
809 }
810 CoolPropDbl d4alphar_dxi_dDelta2_dTau(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) {
811 if (xN_flag == XN_INDEPENDENT) {
812 return HEOS.components[i].EOS().d3alphar_dDelta2_dTau(HEOS.tau(), HEOS.delta());
813 } else {
814 throw ValueError(format("xN_flag is invalid"));
815 }
816 }
817
818 CoolPropDbl d3alphardxidxjdxk(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, std::size_t k,
819 x_N_dependency_flag xN_flag) {
820 return 0;
821 }
822 CoolPropDbl d3alphar_dxi_dxj_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
823 x_N_dependency_flag xN_flag) {
824 return 0;
825 }
826 CoolPropDbl d3alphar_dxi_dxj_dTau(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
827 x_N_dependency_flag xN_flag) {
828 return 0;
829 }
830 CoolPropDbl d4alphar_dxi_dxj_dDelta2(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
831 x_N_dependency_flag xN_flag) {
832 return 0;
833 }
834 CoolPropDbl d4alphar_dxi_dxj_dDelta_dTau(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
835 x_N_dependency_flag xN_flag) {
836 return 0;
837 }
838 CoolPropDbl d4alphar_dxi_dxj_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
839 x_N_dependency_flag xN_flag) {
840 return 0;
841 }
842};
843
849class ResidualHelmholtz
850{
851 public:
852 ExcessTerm Excess;
853 CorrespondingStatesTerm CS;
854
855 ResidualHelmholtz() = default;
856 ResidualHelmholtz(const ExcessTerm& E, const CorrespondingStatesTerm& C) : Excess(E), CS(C) {};
857 virtual ~ResidualHelmholtz() = default;
858
859 ResidualHelmholtz copy() {
860 return ResidualHelmholtz(Excess.copy(), CS);
861 }
862 ResidualHelmholtz* copy_ptr() {
863 return new ResidualHelmholtz(Excess.copy(), CS);
864 }
865
866 virtual HelmholtzDerivatives all(HelmholtzEOSMixtureBackend& HEOS, const std::vector<CoolPropDbl>& mole_fractions, double tau, double delta,
867 bool cache_values = false) {
868 HelmholtzDerivatives a = CS.all(HEOS, tau, delta, mole_fractions, cache_values) + Excess.all(tau, delta, mole_fractions, cache_values);
869 a.delta_x_dalphar_ddelta = delta * a.dalphar_ddelta;
870 a.tau_x_dalphar_dtau = tau * a.dalphar_dtau;
871
872 a.delta2_x_d2alphar_ddelta2 = POW2(delta) * a.d2alphar_ddelta2;
873 a.deltatau_x_d2alphar_ddelta_dtau = delta * tau * a.d2alphar_ddelta_dtau;
874 a.tau2_x_d2alphar_dtau2 = POW2(tau) * a.d2alphar_dtau2;
875
876 return a;
877 }
878 virtual CoolPropDbl dalphar_dxi(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
879 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
880 return CS.dalphar_dxi(HEOS, mole_fractions, i, xN_flag) + Excess.dalphar_dxi(mole_fractions, i, xN_flag);
881 }
882 virtual CoolPropDbl d2alphardxidxj(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
883 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
884 return CS.d2alphardxidxj(HEOS, mole_fractions, i, j, xN_flag) + Excess.d2alphardxidxj(mole_fractions, i, j, xN_flag);
885 }
886 virtual CoolPropDbl d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
887 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
888 return CS.d2alphar_dxi_dTau(HEOS, mole_fractions, i, xN_flag) + Excess.d2alphar_dxi_dTau(mole_fractions, i, xN_flag);
889 }
890 virtual CoolPropDbl d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
891 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
892 return CS.d2alphar_dxi_dDelta(HEOS, mole_fractions, i, xN_flag) + Excess.d2alphar_dxi_dDelta(mole_fractions, i, xN_flag);
893 }
894 virtual CoolPropDbl d3alphar_dxi_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
895 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
896 return CS.d3alphar_dxi_dTau2(HEOS, mole_fractions, i, xN_flag) + Excess.d3alphar_dxi_dTau2(mole_fractions, i, xN_flag);
897 }
898 virtual CoolPropDbl d3alphar_dxi_dDelta_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
899 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
900 return CS.d3alphar_dxi_dDelta_dTau(HEOS, mole_fractions, i, xN_flag) + Excess.d3alphar_dxi_dDelta_dTau(mole_fractions, i, xN_flag);
901 }
902 virtual CoolPropDbl d3alphar_dxi_dDelta2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
903 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
904 return CS.d3alphar_dxi_dDelta2(HEOS, mole_fractions, i, xN_flag) + Excess.d3alphar_dxi_dDelta2(mole_fractions, i, xN_flag);
905 }
906 virtual CoolPropDbl d3alphar_dxi_dxj_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
907 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
908 return CS.d3alphar_dxi_dxj_dTau(HEOS, mole_fractions, i, j, xN_flag) + Excess.d3alphar_dxi_dxj_dTau(mole_fractions, i, j, xN_flag);
909 }
910 virtual CoolPropDbl d3alphar_dxi_dxj_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
911 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
912 return CS.d3alphar_dxi_dxj_dDelta(HEOS, mole_fractions, i, j, xN_flag) + Excess.d3alphar_dxi_dxj_dDelta(mole_fractions, i, j, xN_flag);
913 }
914 virtual CoolPropDbl d3alphardxidxjdxk(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
915 x_N_dependency_flag xN_flag) {
916 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
917 return CS.d3alphardxidxjdxk(HEOS, mole_fractions, i, j, k, xN_flag) + Excess.d3alphardxidxjdxk(mole_fractions, i, j, k, xN_flag);
918 }
919
920 virtual CoolPropDbl d4alphar_dxi_dTau3(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
921 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
922 return CS.d4alphar_dxi_dTau3(HEOS, mole_fractions, i, xN_flag) + Excess.d4alphar_dxi_dTau3(mole_fractions, i, xN_flag);
923 }
924 virtual CoolPropDbl d4alphar_dxi_dDelta2_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
925 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
926 return CS.d4alphar_dxi_dDelta2_dTau(HEOS, mole_fractions, i, xN_flag) + Excess.d4alphar_dxi_dDelta2_dTau(mole_fractions, i, xN_flag);
927 }
928 virtual CoolPropDbl d4alphar_dxi_dDelta_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
929 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
930 return CS.d4alphar_dxi_dDelta_dTau2(HEOS, mole_fractions, i, xN_flag) + Excess.d4alphar_dxi_dDelta_dTau2(mole_fractions, i, xN_flag);
931 }
932 virtual CoolPropDbl d4alphar_dxi_dDelta3(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, x_N_dependency_flag xN_flag) {
933 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
934 return CS.d4alphar_dxi_dDelta3(HEOS, mole_fractions, i, xN_flag) + Excess.d4alphar_dxi_dDelta3(mole_fractions, i, xN_flag);
935 }
936 virtual CoolPropDbl d4alphar_dxi_dxj_dTau2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
937 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
938 return CS.d4alphar_dxi_dxj_dTau2(HEOS, mole_fractions, i, j, xN_flag) + Excess.d4alphar_dxi_dxj_dTau2(mole_fractions, i, j, xN_flag);
939 }
940 virtual CoolPropDbl d4alphar_dxi_dxj_dDelta_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
941 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
942 return CS.d4alphar_dxi_dxj_dDelta_dTau(HEOS, mole_fractions, i, j, xN_flag)
943 + Excess.d4alphar_dxi_dxj_dDelta_dTau(mole_fractions, i, j, xN_flag);
944 }
945 virtual CoolPropDbl d4alphar_dxi_dxj_dDelta2(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) {
946 std::vector<CoolPropDbl>& mole_fractions = HEOS.get_mole_fractions_ref();
947 return CS.d4alphar_dxi_dxj_dDelta2(HEOS, mole_fractions, i, j, xN_flag) + Excess.d4alphar_dxi_dxj_dDelta2(mole_fractions, i, j, xN_flag);
948 }
949 virtual CoolPropDbl d4alphar_dxi_dxj_dxk_dDelta(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
950 x_N_dependency_flag xN_flag) {
951 return 0;
952 }
953 virtual CoolPropDbl d4alphar_dxi_dxj_dxk_dTau(HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
954 x_N_dependency_flag xN_flag) {
955 return 0;
956 }
957};
958
959} /* namespace CoolProp */
960#endif /* HELMHOLTZEOSMIXTUREBACKEND_H_ */
Generated by doxygen 1.9.4