CoolProp 7.2.0
An open-source fluid property and humid air property database
AbstractState.cpp
Go to the documentation of this file.
1/*
2 * AbstractState.cpp
3 *
4 * Created on: 21 Dec 2013
5 * Author: jowr
6 */
7
8#ifndef _CRT_SECURE_NO_WARNINGS
9#define _CRT_SECURE_NO_WARNINGS
10#endif
11
12#include <cstdlib>
13#include <cmath>
14#include "AbstractState.h"
15#include "DataStructures.h"
16#include "Backends/IF97/IF97Backend.h"
17#include "Backends/Cubics/CubicBackend.h"
18#include "Backends/Cubics/VTPRBackend.h"
19#include "Backends/Incompressible/IncompressibleBackend.h"
20#include "Backends/PCSAFT/PCSAFTBackend.h"
21
22#if !defined(NO_TABULAR_BACKENDS)
23#include "Backends/Tabular/TTSEBackend.h"
24#include "Backends/Tabular/BicubicBackend.h"
25#endif
26
27namespace CoolProp {
28
31
32class BackendLibrary
33{
34 private:
35 std::map<backend_families, shared_ptr<AbstractStateGenerator>> backends;
36
37 public:
38 void add_backend(const backend_families& bg, const shared_ptr<AbstractStateGenerator>& asg) {
39 backends[bg] = asg;
40 };
41 void get_generator_iterators(const backend_families& bg,
42 std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator& generator,
43 std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator& end) {
44 generator = backends.find(bg);
45 end = backends.end();
46 };
47 std::size_t size() {
48 return backends.size();
49 };
50};
51inline BackendLibrary& get_backend_library() {
52 static BackendLibrary the_library;
53 return the_library;
54}
55
56void register_backend(const backend_families& bf, shared_ptr<AbstractStateGenerator> gen) {
57 get_backend_library().add_backend(bf, gen);
58};
59
60class IF97BackendGenerator : public AbstractStateGenerator
61{
62 public:
63 AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
64 if (fluid_names.size() == 1) { // Check that fluid_names[0] has only one component
65 const std::string str = fluid_names[0]; // Check that the fluid name is an alias for "Water"
66 if ((upper(str) == "WATER") || (upper(str) == "H2O")) {
67 return new IF97Backend();
68 } else {
69 throw ValueError(format("The IF97 backend returns Water props only; fluid name [%s] not allowed", fluid_names[0].c_str()));
70 }
71 } else {
72 throw ValueError(format("The IF97 backend does not support mixtures, only Water"));
73 };
74 };
75};
76// This static initialization will cause the generator to register
77static GeneratorInitializer<IF97BackendGenerator> if97_gen(IF97_BACKEND_FAMILY);
78class SRKGenerator : public AbstractStateGenerator
79{
80 public:
81 AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
82 return new SRKBackend(fluid_names, get_config_double(R_U_CODATA));
83 };
84};
85static GeneratorInitializer<SRKGenerator> srk_gen(CoolProp::SRK_BACKEND_FAMILY);
86class PRGenerator : public AbstractStateGenerator
87{
88 public:
89 AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
90 return new PengRobinsonBackend(fluid_names, get_config_double(R_U_CODATA));
91 };
92};
93static GeneratorInitializer<PRGenerator> pr_gen(CoolProp::PR_BACKEND_FAMILY);
94class IncompressibleBackendGenerator : public AbstractStateGenerator
95{
96 public:
97 AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
98 if (fluid_names.size() != 1) {
99 throw ValueError(format("For INCOMP backend, name vector must be one element long"));
100 }
101 return new IncompressibleBackend(fluid_names[0]);
102 };
103};
104// This static initialization will cause the generator to register
105static GeneratorInitializer<IncompressibleBackendGenerator> incomp_gen(INCOMP_BACKEND_FAMILY);
106class VTPRGenerator : public CoolProp::AbstractStateGenerator
107{
108 public:
109 CoolProp::AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
110 return new CoolProp::VTPRBackend(fluid_names, CoolProp::get_config_double(R_U_CODATA));
111 };
112};
113// This static initialization will cause the generator to register
114static CoolProp::GeneratorInitializer<VTPRGenerator> vtpr_gen(CoolProp::VTPR_BACKEND_FAMILY);
115
116class PCSAFTGenerator : public CoolProp::AbstractStateGenerator
117{
118 public:
119 CoolProp::AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
120 return new CoolProp::PCSAFTBackend(fluid_names);
121 };
122};
123// This static initialization will cause the generator to register
124static CoolProp::GeneratorInitializer<PCSAFTGenerator> pcsaft_gen(CoolProp::PCSAFT_BACKEND_FAMILY);
125
126AbstractState* AbstractState::factory(const std::string& backend, const std::vector<std::string>& fluid_names) {
127 if (get_debug_level() > 0) {
128 std::cout << "AbstractState::factory(" << backend << "," << stringvec_to_string(fluid_names) << ")" << std::endl;
129 }
130
131 backend_families f1;
132 std::string f2;
133 extract_backend_families_string(backend, f1, f2);
134
135 std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator gen, end;
136 get_backend_library().get_generator_iterators(f1, gen, end);
137
138 if (get_debug_level() > 0) {
139 std::cout << "AbstractState::factory backend_library size: " << get_backend_library().size() << std::endl;
140 }
141
142 if (gen != end) {
143 // One of the registered backends was able to match the given backend family
144 return gen->second->get_AbstractState(fluid_names);
145 }
146#if !defined(NO_TABULAR_BACKENDS)
147 else if (f1 == TTSE_BACKEND_FAMILY) {
148 // Will throw if there is a problem with this backend
149 const shared_ptr<AbstractState> AS(factory(f2, fluid_names));
150 return new TTSEBackend(AS);
151 } else if (f1 == BICUBIC_BACKEND_FAMILY) {
152 // Will throw if there is a problem with this backend
153 const shared_ptr<AbstractState> AS(factory(f2, fluid_names));
154 return new BicubicBackend(AS);
155 }
156#endif
157 else if (backend == "?" || backend.empty()) {
158 const std::size_t idel = fluid_names[0].find("::");
159 // Backend has not been specified, and we have to figure out what the backend is by parsing the string
160 if (idel == std::string::npos) // No '::' found, no backend specified, try HEOS, otherwise a failure
161 {
162 // Figure out what backend to use
163 return factory("HEOS", fluid_names);
164 } else {
165 // Split string at the '::' into two std::string, call again
166 return factory(std::string(fluid_names[0].begin(), fluid_names[0].begin() + idel),
167 std::string(fluid_names[0].begin() + (idel + 2), fluid_names[0].end()));
168 }
169 } else {
170 throw ValueError(format("Invalid backend name [%s] to factory function", backend.c_str()));
171 }
172}
173std::vector<std::string> AbstractState::fluid_names() {
174 return calc_fluid_names();
175}
176bool AbstractState::clear_comp_change() {
177 // Reset all instances of CachedElement and overwrite
178 // the internal double values with -_HUGE
179 this->_molar_mass.clear();
180 this->_critical.fill(_HUGE);
181 this->_reducing.fill(_HUGE);
182 calc_reducing_state();
183 this->_gas_constant = calc_gas_constant();
184 return true;
185}
186bool AbstractState::clear() {
187 // Reset all instances of CachedElement
188 cache.clear();
189
190 this->_critical.fill(_HUGE);
191
193 this->_rhomolar = -_HUGE;
194 this->_T = -_HUGE;
195 this->_p = -_HUGE;
196 this->_Q = -_HUGE;
197
198 return true;
199}
200void AbstractState::mass_to_molar_inputs(CoolProp::input_pairs& input_pair, CoolPropDbl& value1, CoolPropDbl& value2) {
201 // Check if a mass based input, convert it to molar units
202
203 switch (input_pair) {
204 case DmassT_INPUTS:
205 //case HmassT_INPUTS: ///< Enthalpy in J/kg, Temperature in K (NOT CURRENTLY IMPLEMENTED)
206 case SmassT_INPUTS:
207 //case TUmass_INPUTS: ///< Temperature in K, Internal energy in J/kg (NOT CURRENTLY IMPLEMENTED)
208 case DmassP_INPUTS:
209 case DmassQ_INPUTS:
210 case HmassP_INPUTS:
211 case PSmass_INPUTS:
212 case PUmass_INPUTS:
213 case HmassSmass_INPUTS:
214 case SmassUmass_INPUTS:
215 case DmassHmass_INPUTS:
216 case DmassSmass_INPUTS:
217 case DmassUmass_INPUTS:
218 {
219 // Set the cache value for the molar mass if it hasn't been set yet
220 molar_mass();
221
222 // Molar mass (just for compactness of the following switch)
223 const auto mm = static_cast<CoolPropDbl>(_molar_mass);
224
225 switch (input_pair) {
226 case DmassT_INPUTS:
227 input_pair = DmolarT_INPUTS;
228 value1 /= mm;
229 break;
230 //case HmassT_INPUTS: input_pair = HmolarT_INPUTS; value1 *= mm; break; (NOT CURRENTLY IMPLEMENTED)
231 case SmassT_INPUTS:
232 input_pair = SmolarT_INPUTS;
233 value1 *= mm;
234 break;
235 //case TUmass_INPUTS: input_pair = TUmolar_INPUTS; value2 *= mm; break; (NOT CURRENTLY IMPLEMENTED)
236 case DmassP_INPUTS:
237 input_pair = DmolarP_INPUTS;
238 value1 /= mm;
239 break;
240 case DmassQ_INPUTS:
241 input_pair = DmolarQ_INPUTS;
242 value1 /= mm;
243 break;
244 case HmassP_INPUTS:
245 input_pair = HmolarP_INPUTS;
246 value1 *= mm;
247 break;
248 case PSmass_INPUTS:
249 input_pair = PSmolar_INPUTS;
250 value2 *= mm;
251 break;
252 case PUmass_INPUTS:
253 input_pair = PUmolar_INPUTS;
254 value2 *= mm;
255 break;
256 case HmassSmass_INPUTS:
257 input_pair = HmolarSmolar_INPUTS;
258 value1 *= mm;
259 value2 *= mm;
260 break;
261 case SmassUmass_INPUTS:
262 input_pair = SmolarUmolar_INPUTS;
263 value1 *= mm;
264 value2 *= mm;
265 break;
266 case DmassHmass_INPUTS:
267 input_pair = DmolarHmolar_INPUTS;
268 value1 /= mm;
269 value2 *= mm;
270 break;
271 case DmassSmass_INPUTS:
272 input_pair = DmolarSmolar_INPUTS;
273 value1 /= mm;
274 value2 *= mm;
275 break;
276 case DmassUmass_INPUTS:
277 input_pair = DmolarUmolar_INPUTS;
278 value1 /= mm;
279 value2 *= mm;
280 break;
281 default:
282 break;
283 }
284 break;
285 }
286 default:
287 return;
288 }
289}
290double AbstractState::trivial_keyed_output(parameters key) {
291 if (get_debug_level() >= 50)
292 std::cout << format("AbstractState: trivial_keyed_output called for %s ", get_parameter_information(key, "short").c_str()) << std::endl;
293 switch (key) {
294 case imolar_mass:
295 return molar_mass();
296 case iacentric_factor:
297 return acentric_factor();
298 case igas_constant:
299 return gas_constant();
300 case iT_min:
301 return Tmin();
302 case iT_triple:
303 return Ttriple();
304 case iT_max:
305 return Tmax();
306 case iP_max:
307 return pmax();
308 case iP_min:
309 case iP_triple:
310 return this->p_triple();
311 case iT_reducing:
312 return calc_T_reducing();
313 case irhomolar_reducing:
314 return calc_rhomolar_reducing();
315 case iP_reducing:
316 return calc_p_reducing();
317 case iP_critical:
318 return this->p_critical();
319 case iT_critical:
320 return this->T_critical();
321 case irhomolar_critical:
322 return this->rhomolar_critical();
323 case irhomass_critical:
324 return this->rhomass_critical();
325 case iODP:
326 return this->calc_ODP();
327 case iGWP100:
328 return this->calc_GWP100();
329 case iGWP20:
330 return this->calc_GWP20();
331 case iGWP500:
332 return this->calc_GWP500();
333 case ifraction_min:
334 return this->calc_fraction_min();
335 case ifraction_max:
336 return this->calc_fraction_max();
337 case iT_freeze:
338 return this->calc_T_freeze();
339 case iFH:
340 return this->calc_flame_hazard();
341 case iHH:
342 return this->calc_health_hazard();
343 case iPH:
344 return this->calc_physical_hazard();
345 case idipole_moment:
346 return this->calc_dipole_moment();
347 default:
348 throw ValueError(
349 format("This input [%d: \"%s\"] is not valid for trivial_keyed_output", key, get_parameter_information(key, "short").c_str()));
350 }
351}
352double AbstractState::keyed_output(parameters key) {
353 if (get_debug_level() >= 50)
354 std::cout << format("AbstractState: keyed_output called for %s ", get_parameter_information(key, "short").c_str()) << std::endl;
355 // Handle trivial inputs
356 if (is_trivial_parameter(key)) {
357 return trivial_keyed_output(key);
358 }
359 switch (key) {
360 case iQ:
361 return Q();
362 case iT:
363 return T();
364 case iP:
365 return p();
366 case iDmolar:
367 return rhomolar();
368 case iDmass:
369 return rhomass();
370 case iHmolar:
371 return hmolar();
372 case iHmolar_residual:
373 return hmolar_residual();
374 case iHmolar_idealgas:
375 return hmolar_idealgas();
376 case iHmass:
377 return hmass();
378 case iHmass_idealgas:
379 return hmass_idealgas();
380 case iSmolar:
381 return smolar();
382 case iSmolar_residual:
383 return smolar_residual();
384 case iSmolar_idealgas:
385 return smolar_idealgas();
386 case iSmass:
387 return smass();
388 case iSmass_idealgas:
389 return smass_idealgas();
390 case iUmolar:
391 return umolar();
392 case iUmolar_idealgas:
393 return umolar_idealgas();
394 case iUmass:
395 return umass();
396 case iUmass_idealgas:
397 return umass_idealgas();
398 case iGmolar:
399 return gibbsmolar();
400 case iGmolar_residual:
401 return gibbsmolar_residual();
402 case iGmass:
403 return gibbsmass();
404 case iHelmholtzmolar:
405 return helmholtzmolar();
406 case iHelmholtzmass:
407 return helmholtzmass();
408 case iCvmolar:
409 return cvmolar();
410 case iCvmass:
411 return cvmass();
412 case iCpmolar:
413 return cpmolar();
414 case iCp0molar:
415 return cp0molar();
416 case iCpmass:
417 return cpmass();
418 case iCp0mass:
419 return cp0mass();
420 case imolar_mass:
421 return molar_mass();
422 case iT_reducing:
423 return get_reducing_state().T;
424 case irhomolar_reducing:
425 return get_reducing_state().rhomolar;
426 case ispeed_sound:
427 return speed_sound();
428 case ialphar:
429 return alphar();
430 case ialpha0:
431 return alpha0();
432 case idalpha0_ddelta_consttau:
433 return dalpha0_dDelta();
434 case id2alpha0_ddelta2_consttau:
435 return d2alpha0_dDelta2();
436 case id3alpha0_ddelta3_consttau:
437 return d3alpha0_dDelta3();
438 case idalpha0_dtau_constdelta:
439 return dalpha0_dTau();
440 case idalphar_ddelta_consttau:
441 return dalphar_dDelta();
442 case idalphar_dtau_constdelta:
443 return dalphar_dTau();
444 case iBvirial:
445 return Bvirial();
446 case idBvirial_dT:
447 return dBvirial_dT();
448 case iCvirial:
449 return Cvirial();
450 case idCvirial_dT:
451 return dCvirial_dT();
452 case iisothermal_compressibility:
453 return isothermal_compressibility();
454 case iisobaric_expansion_coefficient:
455 return isobaric_expansion_coefficient();
456 case iisentropic_expansion_coefficient:
457 return isentropic_expansion_coefficient();
458 case iviscosity:
459 return viscosity();
460 case iconductivity:
461 return conductivity();
462 case iPrandtl:
463 return Prandtl();
464 case isurface_tension:
465 return surface_tension();
466 case iPhase:
467 return phase();
468 case iZ:
469 return compressibility_factor();
470 case iPIP:
471 return PIP();
472 case ifundamental_derivative_of_gas_dynamics:
473 return fundamental_derivative_of_gas_dynamics();
474 case iTau:
475 return _reducing.T / _T;
476 case iDelta:
477 return _rhomolar / _reducing.rhomolar;
478 default:
479 throw ValueError(format("This input [%d: \"%s\"] is not valid for keyed_output", key, get_parameter_information(key, "short").c_str()));
480 }
481}
482
483double AbstractState::tau() {
484 if (!_tau) _tau = calc_reciprocal_reduced_temperature();
485 return _tau;
486}
487double AbstractState::delta() {
488 if (!_delta) _delta = calc_reduced_density();
489 return _delta;
490}
491double AbstractState::Tmin() {
492 return calc_Tmin();
493}
494double AbstractState::Tmax() {
495 return calc_Tmax();
496}
497double AbstractState::Ttriple() {
498 return calc_Ttriple();
499}
500double AbstractState::pmax() {
501 return calc_pmax();
502}
503double AbstractState::T_critical() {
504 return calc_T_critical();
505}
506double AbstractState::T_reducing() {
507 if (!ValidNumber(_reducing.T)) {
508 calc_reducing_state();
509 }
510 return _reducing.T;
511}
512double AbstractState::p_critical() {
513 return calc_p_critical();
514}
515double AbstractState::p_triple() {
516 return calc_p_triple();
517}
518double AbstractState::rhomolar_critical() {
519 return calc_rhomolar_critical();
520}
521double AbstractState::rhomass_critical() {
522 return calc_rhomolar_critical() * molar_mass();
523}
524double AbstractState::rhomolar_reducing() {
525 if (!ValidNumber(_reducing.rhomolar)) {
526 calc_reducing_state();
527 }
528 return _reducing.rhomolar;
529}
530double AbstractState::rhomass_reducing() {
531 return rhomolar_reducing() * molar_mass();
532}
533double AbstractState::hmolar() {
534 if (!_hmolar) _hmolar = calc_hmolar();
535 return _hmolar;
536}
537double AbstractState::hmolar_residual() {
538 if (!_hmolar_residual) _hmolar_residual = calc_hmolar_residual();
539 return _hmolar_residual;
540}
541double AbstractState::hmolar_idealgas() {
542 return gas_constant() * T() * (1 + tau() * dalpha0_dTau());
543}
544double AbstractState::hmass_idealgas() {
545 return hmolar_idealgas() / molar_mass();
546}
547double AbstractState::hmolar_excess() {
548 if (!_hmolar_excess) calc_excess_properties();
549 return _hmolar_excess;
550}
551double AbstractState::smolar() {
552 if (!_smolar) _smolar = calc_smolar();
553 return _smolar;
554}
555double AbstractState::smolar_residual() {
556 if (!_smolar_residual) _smolar_residual = calc_smolar_residual();
557 return _smolar_residual;
558}
559double AbstractState::smolar_idealgas() {
560 return gas_constant() * (tau() * dalpha0_dTau() - alpha0());
561}
562double AbstractState::smass_idealgas() {
563 return smolar_idealgas() / molar_mass();
564}
565double AbstractState::neff() {
566 const double tau = calc_T_reducing() / _T;
567 const double delta = _rhomolar / calc_rhomolar_reducing();
568 const double Ar01 = delta * dalphar_dDelta();
569 const double Ar11 = tau * delta * d2alphar_dDelta_dTau();
570 const double Ar20 = tau * tau * d2alphar_dTau2();
571 return -3.0 * (Ar01 - Ar11) / Ar20;
572}
573double AbstractState::smolar_excess() {
574 if (!_smolar_excess) calc_excess_properties();
575 return _smolar_excess;
576}
577double AbstractState::umolar() {
578 if (!_umolar) _umolar = calc_umolar();
579 return _umolar;
580}
581double AbstractState::umolar_excess() {
582 if (!_umolar_excess) calc_excess_properties();
583 return _umolar_excess;
584}
585double AbstractState::umolar_idealgas() {
586 return gas_constant() * T() * (tau() * dalpha0_dTau());
587}
588double AbstractState::umass_idealgas() {
589 return umolar_idealgas() / molar_mass();
590}
591double AbstractState::gibbsmolar() {
592 if (!_gibbsmolar) _gibbsmolar = calc_gibbsmolar();
593 return _gibbsmolar;
594}
595double AbstractState::gibbsmolar_residual() {
596 if (!_gibbsmolar_residual) _gibbsmolar_residual = calc_gibbsmolar_residual();
597 return _gibbsmolar_residual;
598}
599double AbstractState::gibbsmolar_excess() {
600 if (!_gibbsmolar_excess) calc_excess_properties();
601 return _gibbsmolar_excess;
602}
603double AbstractState::helmholtzmolar() {
604 if (!_helmholtzmolar) _helmholtzmolar = calc_helmholtzmolar();
605 return _helmholtzmolar;
606}
607double AbstractState::helmholtzmolar_excess() {
608 if (!_helmholtzmolar_excess) calc_excess_properties();
609 return _helmholtzmolar_excess;
610}
611double AbstractState::volumemolar_excess() {
612 if (!_volumemolar_excess) calc_excess_properties();
613 return _volumemolar_excess;
614}
615double AbstractState::cpmolar() {
616 if (!_cpmolar) _cpmolar = calc_cpmolar();
617 return _cpmolar;
618}
619double AbstractState::cp0molar() {
620 return calc_cpmolar_idealgas();
621}
622double AbstractState::cvmolar() {
623 if (!_cvmolar) _cvmolar = calc_cvmolar();
624 return _cvmolar;
625}
626double AbstractState::speed_sound() {
627 if (!_speed_sound) _speed_sound = calc_speed_sound();
628 return _speed_sound;
629}
630double AbstractState::viscosity() {
631 if (!_viscosity) _viscosity = calc_viscosity();
632 return _viscosity;
633}
634double AbstractState::conductivity() {
635 if (!_conductivity) _conductivity = calc_conductivity();
636 return _conductivity;
637}
638double AbstractState::melting_line(int param, int given, double value) {
639 return calc_melting_line(param, given, value);
640}
641double AbstractState::acentric_factor() {
642 return calc_acentric_factor();
643}
644double AbstractState::saturation_ancillary(parameters param, int Q, parameters given, double value) {
645 return calc_saturation_ancillary(param, Q, given, value);
646}
647double AbstractState::surface_tension() {
648 if (!_surface_tension) _surface_tension = calc_surface_tension();
649 return _surface_tension;
650}
651double AbstractState::molar_mass() {
652 if (!_molar_mass) _molar_mass = calc_molar_mass();
653 return _molar_mass;
654}
655double AbstractState::gas_constant() {
656 if (!_gas_constant) _gas_constant = calc_gas_constant();
657 return _gas_constant;
658}
659double AbstractState::fugacity_coefficient(std::size_t i) {
660 // TODO: Cache the fug. coeff for each component
661 return calc_fugacity_coefficient(i);
662}
663std::vector<double> AbstractState::fugacity_coefficients() {
664 // TODO: Cache the fug. coeff for each component
665 return calc_fugacity_coefficients();
666}
667double AbstractState::fugacity(std::size_t i) {
668 // TODO: Cache the fug. coeff for each component
669 return calc_fugacity(i);
670}
671double AbstractState::chemical_potential(std::size_t i) {
672 // TODO: Cache the chemical potential for each component
673 return calc_chemical_potential(i);
674}
675void AbstractState::build_phase_envelope(const std::string& type) {
676 calc_phase_envelope(type);
677}
678double AbstractState::isothermal_compressibility() {
679 return 1.0 / _rhomolar * first_partial_deriv(iDmolar, iP, iT);
680}
681double AbstractState::isobaric_expansion_coefficient() {
682 return -1.0 / _rhomolar * first_partial_deriv(iDmolar, iT, iP);
683}
684double AbstractState::isentropic_expansion_coefficient() {
685 return _rhomolar / _p * first_partial_deriv(iP, iDmolar, iSmolar);
686}
687double AbstractState::Bvirial() {
688 return calc_Bvirial();
689}
690double AbstractState::Cvirial() {
691 return calc_Cvirial();
692}
693double AbstractState::dBvirial_dT() {
694 return calc_dBvirial_dT();
695}
696double AbstractState::dCvirial_dT() {
697 return calc_dCvirial_dT();
698}
699double AbstractState::compressibility_factor() {
700 return calc_compressibility_factor();
701}
702
703double AbstractState::fundamental_derivative_of_gas_dynamics() {
704 // See Colonna, FPE, 2010, Eq. 1
705 return 1 + this->second_partial_deriv(iP, iDmass, iSmolar, iDmass, iSmolar) * this->rhomass() / (2 * powInt(speed_sound(), 2));
706};
707
708// Get the derivatives of the parameters in the partial derivative with respect to T and rho
709void get_dT_drho(AbstractState& AS, parameters index, CoolPropDbl& dT, CoolPropDbl& drho) {
710 const CoolPropDbl T = AS.T(), rho = AS.rhomolar(), rhor = AS.rhomolar_reducing(), Tr = AS.T_reducing(), dT_dtau = -pow(T, 2) / Tr,
711 R = AS.gas_constant(), delta = rho / rhor, tau = Tr / T;
712
713 switch (index) {
714 case iT:
715 dT = 1;
716 drho = 0;
717 break;
718 case iDmolar:
719 dT = 0;
720 drho = 1;
721 break;
722 case iDmass:
723 dT = 0;
724 drho = AS.molar_mass();
725 break;
726 case iP: {
727 // dp/drho|T
728 drho = R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2());
729 // dp/dT|rho
730 dT = rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau());
731 break;
732 }
733 case iHmass:
734 case iHmolar: {
735 // dh/dT|rho
736 dT = R
737 * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
738 + (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()));
739 // dh/drhomolar|T
740 drho = T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2());
741 if (index == iHmass) {
742 // dhmolar/drhomolar|T * dhmass/dhmolar where dhmass/dhmolar = 1/mole mass
743 drho /= AS.molar_mass();
744 dT /= AS.molar_mass();
745 }
746 break;
747 }
748 case iSmass:
749 case iSmolar: {
750 // ds/dT|rho
751 dT = R / T * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
752 // ds/drho|T
753 drho = R / rho * (-(1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()));
754 if (index == iSmass) {
755 // ds/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
756 drho /= AS.molar_mass();
757 dT /= AS.molar_mass();
758 }
759 break;
760 }
761 case iUmass:
762 case iUmolar: {
763 // du/dT|rho
764 dT = R * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
765 // du/drho|T
766 drho = AS.T() * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau());
767 if (index == iUmass) {
768 // du/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
769 drho /= AS.molar_mass();
770 dT /= AS.molar_mass();
771 }
772 break;
773 }
774 case iGmass:
775 case iGmolar: {
776 // dg/dT|rho
777 const double dTau_dT = 1 / dT_dtau;
778 dT = R * AS.T() * (AS.dalpha0_dTau() + AS.dalphar_dTau() + AS.delta() * AS.d2alphar_dDelta_dTau()) * dTau_dT
779 + R * (1 + AS.alpha0() + AS.alphar() + AS.delta() * AS.dalphar_dDelta());
780 // dg/drho|T
781 const double dDelta_drho = 1 / rhor;
782 drho = AS.T() * R * (AS.dalpha0_dDelta() + AS.dalphar_dDelta() + AS.delta() * AS.d2alphar_dDelta2() + AS.dalphar_dDelta()) * dDelta_drho;
783 if (index == iGmass) {
784 // dg/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
785 drho /= AS.molar_mass();
786 dT /= AS.molar_mass();
787 }
788 break;
789 }
790 case iTau:
791 dT = 1 / dT_dtau;
792 drho = 0;
793 break;
794 case iDelta:
795 dT = 0;
796 drho = 1 / rhor;
797 break;
798 case iCvmolar:
799 case iCvmass: {
800 // use the second order derivative of internal energy
801 // make it cleaner by using the function get_dT_drho_second_derivatives directly?
802 // dcvdT|rho = d2u/dT2|rho
803 dT = R / T * pow(tau, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
804 // dcvdrho|T = d2u/dT/drho
805 drho = R / rho * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
806 if (index == iCvmass) {
807 drho /= AS.molar_mass();
808 dT /= AS.molar_mass();
809 }
810 break;
811 }
812 case iCpmolar:
813 case iCpmass: {
814 // dcp/dT|rho = d2h/dT2 + dh/drho * dP/dT * d2P/drhodT / ( dp/drho )^2 - ( d2h/dTdrho * dP/dT + dh/drho * d2P/dT2 ) / ( dP/drho )
815 dT = R / T * pow(tau, 2)
816 * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
817 + delta * AS.d3alphar_dDelta_dTau2());
818 dT += (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
819 * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
820 * (R
821 * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2() - 2 * delta * tau * AS.d2alphar_dDelta_dTau()
822 - tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau()))
823 / pow(R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()), 2);
824 dT -= ((R / rho * delta
825 * (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
826 - tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau()))
827 * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
828 + (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
829 * (rho * R / T * (pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2())))
830 / (R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()));
831 // dcpdrho|T = d2h/dTdrho + dh/drho * dP/dT * d2P/drho2 / ( dp/drho )^2 - ( d2h/drho2 * dP/dT + dh/drho * d2P/dTdrho ) / ( dP/drho )
832 drho = R / rho * delta
833 * (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
834 - tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau()); //d2h/dTdrho
835 drho +=
836 (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
837 * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
838 * (T * R / rho * (2 * delta * AS.dalphar_dDelta() + 4 * pow(delta, 2) * AS.d2alphar_dDelta2() + pow(delta, 3) * AS.d3alphar_dDelta3()))
839 / pow(R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()), 2);
840 drho -= ((R * T * pow(delta / rho, 2) * (tau * AS.d3alphar_dDelta2_dTau() + 2 * AS.d2alphar_dDelta2() + delta * AS.d3alphar_dDelta3()))
841 * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
842 + (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
843 * (R
844 * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()
845 - 2 * delta * tau * AS.d2alphar_dDelta_dTau() - tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau())))
846 / (R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()));
847 if (index == iCpmass) {
848 drho /= AS.molar_mass();
849 dT /= AS.molar_mass();
850 }
851 break;
852 }
853 case ispeed_sound: {
854 //dwdT
855 const double aa = 1.0 + delta * AS.dalphar_dDelta() - delta * tau * AS.d2alphar_dDelta_dTau();
856 const double bb = pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2());
857 const double daa_dTau = -delta * tau * AS.d3alphar_dDelta_dTau2();
858 const double dbb_dTau =
859 pow(tau, 2) * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2.0 * tau * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2());
860 const double w = AS.speed_sound();
861 dT = 1.0 / 2.0 / w / T
862 * (pow(w, 2)
863 - R * Tr / AS.molar_mass()
864 * (2.0 * delta * AS.d2alphar_dDelta_dTau() + pow(delta, 2) * AS.d3alphar_dDelta2_dTau()
865 - (2 * aa / bb * daa_dTau - pow(aa / bb, 2) * dbb_dTau)));
866 //dwdrho
867 const double daa_dDelta =
868 AS.dalphar_dDelta() + delta * AS.d2alphar_dDelta2() - tau * (AS.d2alphar_dDelta_dTau() + delta * AS.d3alphar_dDelta2_dTau());
869 const double dbb_dDelta = pow(tau, 2) * (AS.d3alpha0_dDelta_dTau2() + AS.d3alphar_dDelta_dTau2());
870 drho = R * T / 2.0 / AS.molar_mass() / w / rhor
871 * (2.0 * (AS.dalphar_dDelta() + delta * AS.d2alphar_dDelta2())
872 + (2.0 * delta * AS.d2alphar_dDelta2() + pow(delta, 2) * AS.d3alphar_dDelta3())
873 - (2 * aa / bb * daa_dDelta - pow(aa / bb, 2) * dbb_dDelta));
874 break;
875 }
876 default:
877 throw ValueError(format("input to get_dT_drho[%s] is invalid", get_parameter_information(index, "short").c_str()));
878 }
879}
880void get_dT_drho_second_derivatives(AbstractState& AS, int index, CoolPropDbl& dT2, CoolPropDbl& drho_dT, CoolPropDbl& drho2) {
881 const CoolPropDbl T = AS.T(), rho = AS.rhomolar(), rhor = AS.rhomolar_reducing(), Tr = AS.T_reducing(), R = AS.gas_constant(), delta = rho / rhor,
882 tau = Tr / T;
883
884 // Here we use T and rho as independent variables since derivations are already done by Thorade, 2013,
885 // Partial derivatives of thermodynamic state propertiesfor dynamic simulation, DOI 10.1007/s12665-013-2394-z
886
887 switch (index) {
888 case iT:
889 case iDmass:
890 case iDmolar:
891 dT2 = 0; // d2rhomolar_dtau2
892 drho2 = 0;
893 drho_dT = 0;
894 break;
895 case iTau:
896 dT2 = 2 * Tr / pow(T, 3);
897 drho_dT = 0;
898 drho2 = 0;
899 break;
900 case iDelta:
901 dT2 = 0;
902 drho_dT = 0;
903 drho2 = 0;
904 break;
905 case iP: {
906 drho2 =
907 T * R / rho * (2 * delta * AS.dalphar_dDelta() + 4 * pow(delta, 2) * AS.d2alphar_dDelta2() + pow(delta, 3) * AS.d3alphar_dDelta3());
908 dT2 = rho * R / T * (pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
909 drho_dT = R
910 * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2() - 2 * delta * tau * AS.d2alphar_dDelta_dTau()
911 - tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau());
912 break;
913 }
914 case iHmass:
915 case iHmolar: {
916 // d2h/drho2|T
917 drho2 = R * T * pow(delta / rho, 2) * (tau * AS.d3alphar_dDelta2_dTau() + 2 * AS.d2alphar_dDelta2() + delta * AS.d3alphar_dDelta3());
918 // d2h/dT2|rho
919 dT2 = R / T * pow(tau, 2)
920 * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
921 + delta * AS.d3alphar_dDelta_dTau2());
922 // d2h/drho/dT
923 drho_dT = R / rho * delta
924 * (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
925 - tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau());
926 if (index == iHmass) {
927 drho2 /= AS.molar_mass();
928 drho_dT /= AS.molar_mass();
929 dT2 /= AS.molar_mass();
930 }
931 break;
932 }
933 case iSmass:
934 case iSmolar: {
935 // d2s/rho2|T
936 drho2 = R / pow(rho, 2) * (1 - pow(delta, 2) * AS.d2alphar_dDelta2() + tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau());
937 // d2s/dT2|rho
938 dT2 = R * pow(tau / T, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 3 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
939 // d2s/drho/dT
940 drho_dT = R / (T * rho) * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
941 if (index == iSmass) {
942 drho2 /= AS.molar_mass();
943 drho_dT /= AS.molar_mass();
944 dT2 /= AS.molar_mass();
945 }
946 break;
947 }
948 case iUmass:
949 case iUmolar: {
950 // d2u/rho2|T
951 drho2 = R * T * tau * pow(delta / rho, 2) * AS.d3alphar_dDelta2_dTau();
952 // d2u/dT2|rho
953 dT2 = R / T * pow(tau, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
954 // d2u/drho/dT
955 drho_dT = R / rho * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
956 if (index == iUmass) {
957 drho2 /= AS.molar_mass();
958 drho_dT /= AS.molar_mass();
959 dT2 /= AS.molar_mass();
960 }
961 break;
962 }
963 default:
964 throw ValueError(format("input to get_dT_drho_second_derivatives[%s] is invalid", get_parameter_information(index, "short").c_str()));
965 }
966}
967CoolPropDbl AbstractState::calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant) {
968 CoolPropDbl dOf_dT, dOf_drho, dWrt_dT, dWrt_drho, dConstant_dT, dConstant_drho;
969
970 get_dT_drho(*this, Of, dOf_dT, dOf_drho);
971 get_dT_drho(*this, Wrt, dWrt_dT, dWrt_drho);
972 get_dT_drho(*this, Constant, dConstant_dT, dConstant_drho);
973
974 return (dOf_dT * dConstant_drho - dOf_drho * dConstant_dT) / (dWrt_dT * dConstant_drho - dWrt_drho * dConstant_dT);
975}
976CoolPropDbl AbstractState::calc_second_partial_deriv(parameters Of1, parameters Wrt1, parameters Constant1, parameters Wrt2, parameters Constant2) {
977 CoolPropDbl dOf1_dT, dOf1_drho, dWrt1_dT, dWrt1_drho, dConstant1_dT, dConstant1_drho, d2Of1_dT2, d2Of1_drhodT, d2Of1_drho2, d2Wrt1_dT2,
978 d2Wrt1_drhodT, d2Wrt1_drho2, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2, dWrt2_dT, dWrt2_drho, dConstant2_dT, dConstant2_drho, N, D,
979 dNdrho__T, dDdrho__T, dNdT__rho, dDdT__rho, dderiv1_drho, dderiv1_dT, second;
980
981 // First and second partials needed for terms involved in first derivative
982 get_dT_drho(*this, Of1, dOf1_dT, dOf1_drho);
983 get_dT_drho(*this, Wrt1, dWrt1_dT, dWrt1_drho);
984 get_dT_drho(*this, Constant1, dConstant1_dT, dConstant1_drho);
985 get_dT_drho_second_derivatives(*this, Of1, d2Of1_dT2, d2Of1_drhodT, d2Of1_drho2);
986 get_dT_drho_second_derivatives(*this, Wrt1, d2Wrt1_dT2, d2Wrt1_drhodT, d2Wrt1_drho2);
987 get_dT_drho_second_derivatives(*this, Constant1, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2);
988
989 // First derivatives of terms involved in the second derivative
990 get_dT_drho(*this, Wrt2, dWrt2_dT, dWrt2_drho);
991 get_dT_drho(*this, Constant2, dConstant2_dT, dConstant2_drho);
992
993 // Numerator and denominator of first partial derivative term
994 N = dOf1_dT * dConstant1_drho - dOf1_drho * dConstant1_dT;
995 D = dWrt1_dT * dConstant1_drho - dWrt1_drho * dConstant1_dT;
996
997 // Derivatives of the numerator and denominator of the first partial derivative term with respect to rho, T held constant
998 // They are of similar form, with Of1 and Wrt1 swapped
999 dNdrho__T = dOf1_dT * d2Constant1_drho2 + d2Of1_drhodT * dConstant1_drho - dOf1_drho * d2Constant1_drhodT - d2Of1_drho2 * dConstant1_dT;
1000 dDdrho__T = dWrt1_dT * d2Constant1_drho2 + d2Wrt1_drhodT * dConstant1_drho - dWrt1_drho * d2Constant1_drhodT - d2Wrt1_drho2 * dConstant1_dT;
1001
1002 // Derivatives of the numerator and denominator of the first partial derivative term with respect to T, rho held constant
1003 // They are of similar form, with Of1 and Wrt1 swapped
1004 dNdT__rho = dOf1_dT * d2Constant1_drhodT + d2Of1_dT2 * dConstant1_drho - dOf1_drho * d2Constant1_dT2 - d2Of1_drhodT * dConstant1_dT;
1005 dDdT__rho = dWrt1_dT * d2Constant1_drhodT + d2Wrt1_dT2 * dConstant1_drho - dWrt1_drho * d2Constant1_dT2 - d2Wrt1_drhodT * dConstant1_dT;
1006
1007 // First partial of first derivative term with respect to T
1008 dderiv1_drho = (D * dNdrho__T - N * dDdrho__T) / pow(D, 2);
1009
1010 // First partial of first derivative term with respect to rho
1011 dderiv1_dT = (D * dNdT__rho - N * dDdT__rho) / pow(D, 2);
1012
1013 // Complete second derivative
1014 second = (dderiv1_dT * dConstant2_drho - dderiv1_drho * dConstant2_dT) / (dWrt2_dT * dConstant2_drho - dWrt2_drho * dConstant2_dT);
1015
1016 return second;
1017}
1018// // ----------------------------------------
1019// // Smoothing functions for density
1020// // ----------------------------------------
1021// /// A smoothed version of the derivative using a spline curve in the region of x=0 to x=xend
1022// virtual double AbstractState::drhodh_constp_smoothed(double xend);
1023// /// A smoothed version of the derivative using a spline curve in the region of x=0 to x=xend
1024// virtual double AbstractState::drhodp_consth_smoothed(double xend);
1025// /// Density corresponding to the smoothed derivatives in the region of x=0 to x=xend
1026// virtual void AbstractState::rho_smoothed(double xend, double *rho_spline, double *dsplinedh, double *dsplinedp);
1027
1028} /* namespace CoolProp */
1029
1030#ifdef ENABLE_CATCH
1031
1032#include <catch2/catch_all.hpp>
1033
1034TEST_CASE("Check AbstractState", "[AbstractState]") {
1035 SECTION("bad backend") {
1036 CHECK_THROWS(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("DEFINITELY_A_BAD_BACKEND", "Water")));
1037 }
1038 SECTION("good backend - bad fluid") {
1039 CHECK_THROWS(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("HEOS", "DEFINITELY_A_BAD_FLUID")));
1040 }
1041 SECTION("good backend - helmholtz") {
1042 CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("HEOS", "Water")));
1043 }
1044 SECTION("good backend - incomp") {
1045 CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("INCOMP", "DEB")));
1046 }
1047 SECTION("good backend - REFPROP") {
1048 CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("REFPROP", "Water")));
1049 }
1050}
1051
1052TEST_CASE("Check derivatives in first_partial_deriv", "[derivs_in_first_partial_deriv]") {
1053 shared_ptr<CoolProp::AbstractState> Water(CoolProp::AbstractState::factory("HEOS", "Water"));
1054 shared_ptr<CoolProp::AbstractState> WaterplusT(CoolProp::AbstractState::factory("HEOS", "Water"));
1055 shared_ptr<CoolProp::AbstractState> WaterminusT(CoolProp::AbstractState::factory("HEOS", "Water"));
1056 shared_ptr<CoolProp::AbstractState> Waterplusrho(CoolProp::AbstractState::factory("HEOS", "Water"));
1057 shared_ptr<CoolProp::AbstractState> Waterminusrho(CoolProp::AbstractState::factory("HEOS", "Water"));
1058
1059 double dT = 1e-3, drho = 1;
1060 Water->update(CoolProp::PT_INPUTS, 101325, 300);
1061 WaterplusT->update(CoolProp::DmolarT_INPUTS, Water->rhomolar(), 300 + dT);
1062 WaterminusT->update(CoolProp::DmolarT_INPUTS, Water->rhomolar(), 300 - dT);
1063 Waterplusrho->update(CoolProp::DmolarT_INPUTS, Water->rhomolar() + drho, 300);
1064 Waterminusrho->update(CoolProp::DmolarT_INPUTS, Water->rhomolar() - drho, 300);
1065
1066 // Numerical derivatives
1067 CoolPropDbl dP_dT_num = (WaterplusT->p() - WaterminusT->p()) / (2 * dT);
1068 CoolPropDbl dP_drho_num = (Waterplusrho->p() - Waterminusrho->p()) / (2 * drho);
1069
1070 CoolPropDbl dHmolar_dT_num = (WaterplusT->hmolar() - WaterminusT->hmolar()) / (2 * dT);
1071 CoolPropDbl dHmolar_drho_num = (Waterplusrho->hmolar() - Waterminusrho->hmolar()) / (2 * drho);
1072 CoolPropDbl dHmass_dT_num = (WaterplusT->hmass() - WaterminusT->hmass()) / (2 * dT);
1073 CoolPropDbl dHmass_drho_num = (Waterplusrho->hmass() - Waterminusrho->hmass()) / (2 * drho);
1074
1075 CoolPropDbl dSmolar_dT_num = (WaterplusT->smolar() - WaterminusT->smolar()) / (2 * dT);
1076 CoolPropDbl dSmolar_drho_num = (Waterplusrho->smolar() - Waterminusrho->smolar()) / (2 * drho);
1077 CoolPropDbl dSmass_dT_num = (WaterplusT->smass() - WaterminusT->smass()) / (2 * dT);
1078 CoolPropDbl dSmass_drho_num = (Waterplusrho->smass() - Waterminusrho->smass()) / (2 * drho);
1079
1080 CoolPropDbl dUmolar_dT_num = (WaterplusT->umolar() - WaterminusT->umolar()) / (2 * dT);
1081 CoolPropDbl dUmolar_drho_num = (Waterplusrho->umolar() - Waterminusrho->umolar()) / (2 * drho);
1082 CoolPropDbl dUmass_dT_num = (WaterplusT->umass() - WaterminusT->umass()) / (2 * dT);
1083 CoolPropDbl dUmass_drho_num = (Waterplusrho->umass() - Waterminusrho->umass()) / (2 * drho);
1084
1085 CoolPropDbl dGmolar_dT_num = (WaterplusT->gibbsmolar() - WaterminusT->gibbsmolar()) / (2 * dT);
1086 CoolPropDbl dGmolar_drho_num = (Waterplusrho->gibbsmolar() - Waterminusrho->gibbsmolar()) / (2 * drho);
1087 CoolPropDbl dGmass_dT_num = (WaterplusT->gibbsmass() - WaterminusT->gibbsmass()) / (2 * dT);
1088 CoolPropDbl dGmass_drho_num = (Waterplusrho->gibbsmass() - Waterminusrho->gibbsmass()) / (2 * drho);
1089
1090 CoolPropDbl dCvmolar_dT_num = (WaterplusT->cvmolar() - WaterminusT->cvmolar()) / (2 * dT);
1091 CoolPropDbl dCvmolar_drho_num = (Waterplusrho->cvmolar() - Waterminusrho->cvmolar()) / (2 * drho);
1092 CoolPropDbl dCvmass_dT_num = (WaterplusT->cvmass() - WaterminusT->cvmass()) / (2 * dT);
1093 CoolPropDbl dCvmass_drho_num = (Waterplusrho->cvmass() - Waterminusrho->cvmass()) / (2 * drho);
1094
1095 CoolPropDbl dCpmolar_dT_num = (WaterplusT->cpmolar() - WaterminusT->cpmolar()) / (2 * dT);
1096 CoolPropDbl dCpmolar_drho_num = (Waterplusrho->cpmolar() - Waterminusrho->cpmolar()) / (2 * drho);
1097 CoolPropDbl dCpmass_dT_num = (WaterplusT->cpmass() - WaterminusT->cpmass()) / (2 * dT);
1098 CoolPropDbl dCpmass_drho_num = (Waterplusrho->cpmass() - Waterminusrho->cpmass()) / (2 * drho);
1099
1100 CoolPropDbl dspeed_sound_dT_num = (WaterplusT->speed_sound() - WaterminusT->speed_sound()) / (2 * dT);
1101 CoolPropDbl dspeed_sound_drho_num = (Waterplusrho->speed_sound() - Waterminusrho->speed_sound()) / (2 * drho);
1102
1103 // Pressure
1104 CoolPropDbl dP_dT_analyt, dP_drho_analyt;
1105 CoolProp::get_dT_drho(*Water, CoolProp::iP, dP_dT_analyt, dP_drho_analyt);
1106 // Enthalpy
1107 CoolPropDbl dHmolar_dT_analyt, dHmolar_drho_analyt;
1108 CoolProp::get_dT_drho(*Water, CoolProp::iHmolar, dHmolar_dT_analyt, dHmolar_drho_analyt);
1109 CoolPropDbl dHmass_dT_analyt, dHmass_drho_analyt;
1110 CoolProp::get_dT_drho(*Water, CoolProp::iHmass, dHmass_dT_analyt, dHmass_drho_analyt);
1111 // Entropy
1112 CoolPropDbl dSmolar_dT_analyt, dSmolar_drho_analyt;
1113 CoolProp::get_dT_drho(*Water, CoolProp::iSmolar, dSmolar_dT_analyt, dSmolar_drho_analyt);
1114 CoolPropDbl dSmass_dT_analyt, dSmass_drho_analyt;
1115 CoolProp::get_dT_drho(*Water, CoolProp::iSmass, dSmass_dT_analyt, dSmass_drho_analyt);
1116 // Internal energy
1117 CoolPropDbl dUmolar_dT_analyt, dUmolar_drho_analyt;
1118 CoolProp::get_dT_drho(*Water, CoolProp::iUmolar, dUmolar_dT_analyt, dUmolar_drho_analyt);
1119 CoolPropDbl dUmass_dT_analyt, dUmass_drho_analyt;
1120 CoolProp::get_dT_drho(*Water, CoolProp::iUmass, dUmass_dT_analyt, dUmass_drho_analyt);
1121 // Gibbs
1122 CoolPropDbl dGmolar_dT_analyt, dGmolar_drho_analyt;
1123 CoolProp::get_dT_drho(*Water, CoolProp::iGmolar, dGmolar_dT_analyt, dGmolar_drho_analyt);
1124 CoolPropDbl dGmass_dT_analyt, dGmass_drho_analyt;
1125 CoolProp::get_dT_drho(*Water, CoolProp::iGmass, dGmass_dT_analyt, dGmass_drho_analyt);
1126 // Isochoric heat capacity
1127 CoolPropDbl dCvmolar_dT_analyt, dCvmolar_drho_analyt;
1128 CoolProp::get_dT_drho(*Water, CoolProp::iCvmolar, dCvmolar_dT_analyt, dCvmolar_drho_analyt);
1129 CoolPropDbl dCvmass_dT_analyt, dCvmass_drho_analyt;
1130 CoolProp::get_dT_drho(*Water, CoolProp::iCvmass, dCvmass_dT_analyt, dCvmass_drho_analyt);
1131 // Isobaric heat capacity
1132 CoolPropDbl dCpmolar_dT_analyt, dCpmolar_drho_analyt;
1133 CoolProp::get_dT_drho(*Water, CoolProp::iCpmolar, dCpmolar_dT_analyt, dCpmolar_drho_analyt);
1134 CoolPropDbl dCpmass_dT_analyt, dCpmass_drho_analyt;
1135 CoolProp::get_dT_drho(*Water, CoolProp::iCpmass, dCpmass_dT_analyt, dCpmass_drho_analyt);
1136 // Speed of sound
1137 CoolPropDbl dspeed_sound_dT_analyt, dspeed_sound_drho_analyt;
1138 CoolProp::get_dT_drho(*Water, CoolProp::ispeed_sound, dspeed_sound_dT_analyt, dspeed_sound_drho_analyt);
1139
1140 double eps = 1e-3;
1141
1142 CHECK(std::abs(dP_dT_analyt / dP_dT_num - 1) < eps);
1143 CHECK(std::abs(dP_drho_analyt / dP_drho_num - 1) < eps);
1144
1145 CHECK(std::abs(dHmolar_dT_analyt / dHmolar_dT_num - 1) < eps);
1146 CHECK(std::abs(dHmolar_drho_analyt / dHmolar_drho_num - 1) < eps);
1147 CHECK(std::abs(dHmass_dT_analyt / dHmass_dT_num - 1) < eps);
1148 CHECK(std::abs(dHmass_drho_analyt / dHmass_drho_num - 1) < eps);
1149
1150 CHECK(std::abs(dSmolar_dT_analyt / dSmolar_dT_num - 1) < eps);
1151 CHECK(std::abs(dSmolar_drho_analyt / dSmolar_drho_num - 1) < eps);
1152 CHECK(std::abs(dSmass_dT_analyt / dSmass_dT_num - 1) < eps);
1153 CHECK(std::abs(dSmass_drho_analyt / dSmass_drho_num - 1) < eps);
1154
1155 CHECK(std::abs(dUmolar_dT_analyt / dUmolar_dT_num - 1) < eps);
1156 CHECK(std::abs(dUmolar_drho_analyt / dUmolar_drho_num - 1) < eps);
1157 CHECK(std::abs(dUmass_dT_analyt / dUmass_dT_num - 1) < eps);
1158 CHECK(std::abs(dUmass_drho_analyt / dUmass_drho_num - 1) < eps);
1159
1160 CHECK(std::abs(dGmolar_dT_analyt / dGmolar_dT_num - 1) < eps);
1161 CHECK(std::abs(dGmolar_drho_analyt / dGmolar_drho_num - 1) < eps);
1162 CHECK(std::abs(dGmass_dT_analyt / dGmass_dT_num - 1) < eps);
1163 CHECK(std::abs(dGmass_drho_analyt / dGmass_drho_num - 1) < eps);
1164
1165 CHECK(std::abs(dCvmolar_dT_analyt / dCvmolar_dT_num - 1) < eps);
1166 CHECK(std::abs(dCvmolar_drho_analyt / dCvmolar_drho_num - 1) < eps);
1167 CHECK(std::abs(dCvmass_dT_analyt / dCvmass_dT_num - 1) < eps);
1168 CHECK(std::abs(dCvmass_drho_analyt / dCvmass_drho_num - 1) < eps);
1169
1170 CHECK(std::abs(dCpmolar_dT_analyt / dCpmolar_dT_num - 1) < eps);
1171 CHECK(std::abs(dCpmolar_drho_analyt / dCpmolar_drho_num - 1) < eps);
1172 CHECK(std::abs(dCpmass_dT_analyt / dCpmass_dT_num - 1) < eps);
1173 CHECK(std::abs(dCpmass_drho_analyt / dCpmass_drho_num - 1) < eps);
1174
1175 CHECK(std::abs(dspeed_sound_dT_analyt / dspeed_sound_dT_num - 1) < eps);
1176 CHECK(std::abs(dspeed_sound_drho_analyt / dspeed_sound_drho_num - 1) < eps);
1177}
1178
1179#endif
Generated on Sat Nov 1 2025 19:25:21 for CoolProp by doxygen 1.9.4