CoolProp 8.0.0
An open-source fluid property and humid air property database
IncompressibleFluid.cpp
Go to the documentation of this file.
1
2#include "CoolProp/CoolProp.h"
3#include "CoolProp/DataStructures.h"
4#include "CoolProp/fluids/IncompressibleFluid.h"
5#include "IncompressibleLibrary.h"
6#include <cmath>
7#include "CoolProp/numerics/MatrixMath.h"
8#include "CoolProp/numerics/PolyMath.h"
9#include <Eigen/Core>
10
11constexpr double INCOMP_EPSILON = DBL_EPSILON * 100.0;
12constexpr double INCOMP_DELTA = INCOMP_EPSILON * 10.0;
13
14namespace CoolProp {
15
17
22//IncompressibleFluid::IncompressibleFluid();
23
24void IncompressibleFluid::validate() {
25 return;
26 // TODO: Implement validation function
27
28 // u and s have to be of the polynomial type!
29 //throw NotImplementedError("TODO");
30}
31
32bool IncompressibleFluid::is_pure() {
33 if (density.coeffs.cols() == 1) return true;
34 return false;
35}
36
38double IncompressibleFluid::baseExponential(const IncompressibleData& data, double y, double ybase) {
39 Eigen::VectorXd coeffs = makeVector(data.coeffs);
40 size_t r = coeffs.rows(), c = coeffs.cols();
41 if (strict && (r != 3 || c != 1)) {
42 throw ValueError(format("%s (%d): You have to provide a 3,1 matrix of coefficients, not (%d,%d).", __FILE__, __LINE__, r, c));
43 }
44
45 // Guard the function against zero denominators in exp((double)(coeffs[0] / ((y - ybase) + coeffs[1]) - coeffs[2]))
46 auto fnc = [&](double x) { return exp((double)(coeffs[0] / (x)-coeffs[2])); };
47 double x_den = (y - ybase) + coeffs[1];
48 double x_lo = -INCOMP_EPSILON;
49 double x_hi = +INCOMP_EPSILON;
50 if (x_den < x_lo || x_den > x_hi) {
51 return fnc(x_den);
52 }
53 // ... now we know that we are in the danger zone
54 // step away from the zero-crossing
55 x_lo -= INCOMP_DELTA;
56 x_hi += INCOMP_DELTA;
57 const double f_lo = fnc(x_lo);
58 const double f_hi = fnc(x_hi);
59 // Linearize around the zero-crossing
60 return (f_hi - f_lo) / (x_hi - x_lo) * (x_den - x_lo) + f_lo;
61}
62
64double IncompressibleFluid::baseLogexponential(const IncompressibleData& data, double y, double ybase) {
65 Eigen::VectorXd coeffs = makeVector(data.coeffs);
66 size_t r = coeffs.rows(), c = coeffs.cols();
67 if (strict && (r != 3 || c != 1)) {
68 throw ValueError(format("%s (%d): You have to provide a 3,1 matrix of coefficients, not (%d,%d).", __FILE__, __LINE__, r, c));
69 }
70
71 // Guard the function against zero denominators in exp((double)(log((double)(1.0 / ((y - ybase) + coeffs[0]) + 1.0 / ((y - ybase) + coeffs[0]) / ((y - ybase) + coeffs[0]))) * coeffs[1] + coeffs[2]))
72 auto fnc = [&](double x) { return exp((double)(log((double)(1.0 / (x) + 1.0 / (x) / (x))) * coeffs[1] + coeffs[2])); };
73 double x_den = (y - ybase) + coeffs[0];
74 double x_lo = -INCOMP_EPSILON;
75 double x_hi = +INCOMP_EPSILON;
76 if (x_den < x_lo || x_den > x_hi) {
77 return fnc(x_den);
78 }
79 // ... now we know that we are in the danger zone
80 // step away from the zero-crossing
81 x_lo -= INCOMP_DELTA;
82 x_hi += INCOMP_DELTA;
83 const double f_lo = fnc(x_lo);
84 const double f_hi = fnc(x_hi);
85 // Linearize around the zero-crossing
86 return (f_hi - f_lo) / (x_hi - x_lo) * (x_den - x_lo) + f_lo;
87}
88
89double IncompressibleFluid::basePolyOffset(IncompressibleData data, double y, double z) {
90 size_t r = data.coeffs.rows(), c = data.coeffs.cols();
91 double offset = 0.0;
92 double in = 0.0;
93 Eigen::MatrixXd coeffs;
94 if (r > 0 && c > 0) {
95 offset = data.coeffs(0, 0);
96 if (r == 1 && c > 1) { // row vector -> function of z
97 coeffs = Eigen::MatrixXd(data.coeffs.block(0, 1, r, c - 1));
98 in = z;
99 } else if (r > 1 && c == 1) { // column vector -> function of y
100 coeffs = Eigen::MatrixXd(data.coeffs.block(1, 0, r - 1, c));
101 in = y;
102 } else {
103 throw ValueError(format("%s (%d): You have to provide a vector (1D matrix) of coefficients, not (%d,%d).", __FILE__, __LINE__, r, c));
104 }
105 return poly.evaluate(coeffs, in, 0, offset);
106 }
107 throw ValueError(format("%s (%d): You have to provide a vector (1D matrix) of coefficients, not (%d,%d).", __FILE__, __LINE__, r, c));
108}
109
111double IncompressibleFluid::rho(double T, double p, double x) {
112 switch (density.type) {
113 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
114 return poly.evaluate(density.coeffs, T, x, 0, 0, Tbase, xbase);
115 case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
116 return baseExponential(density, T, 0.0);
117 case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
118 return baseLogexponential(density, T, 0.0);
119 case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
120 return exp(poly.evaluate(density.coeffs, T, x, 0, 0, Tbase, xbase));
121 case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
122 return basePolyOffset(density, T, x);
123 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
124 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?", __FILE__,
125 __LINE__, density.type));
126 default:
127 throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.", __FILE__, __LINE__, density.type));
128 }
129}
130
132double IncompressibleFluid::c(double T, double p, double x) {
133 switch (specific_heat.type) {
134 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
135 //throw NotImplementedError("Here you should implement the polynomial.");
136 return poly.evaluate(specific_heat.coeffs, T, x, 0, 0, Tbase, xbase);
137 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
138 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?", __FILE__,
139 __LINE__, specific_heat.type));
140 default:
141 throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for specific heat.", __FILE__, __LINE__,
142 specific_heat.type));
143 }
144}
145
147double IncompressibleFluid::visc(double T, double p, double x) {
148 switch (viscosity.type) {
149 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
150 return poly.evaluate(viscosity.coeffs, T, x, 0, 0, Tbase, xbase);
151 case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
152 return baseExponential(viscosity, T, 0.0);
153 case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
154 return baseLogexponential(viscosity, T, 0.0);
155 case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
156 return exp(poly.evaluate(viscosity.coeffs, T, x, 0, 0, Tbase, xbase));
157 case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
158 return basePolyOffset(viscosity, T, x);
159 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
160 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?", __FILE__,
161 __LINE__, viscosity.type));
162 default:
163 throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.", __FILE__, __LINE__, viscosity.type));
164 }
165}
167double IncompressibleFluid::cond(double T, double p, double x) {
168 switch (conductivity.type) {
169 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
170 return poly.evaluate(conductivity.coeffs, T, x, 0, 0, Tbase, xbase);
171 case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
172 return baseExponential(conductivity, T, 0.0);
173 case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
174 return baseLogexponential(conductivity, T, 0.0);
175 case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
176 return exp(poly.evaluate(conductivity.coeffs, T, x, 0, 0, Tbase, xbase));
177 case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
178 return basePolyOffset(conductivity, T, x);
179 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
180 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?", __FILE__,
181 __LINE__, conductivity.type));
182 default:
183 throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.", __FILE__, __LINE__, conductivity.type));
184 }
185}
187double IncompressibleFluid::psat(double T, double x) {
188 // Below TminPsat the polynomial fit is not valid. Returning 0.0
189 // silently here surfaced as "Psat = 0 for the entire valid T range"
190 // for fluids whose TminPsat is at or above Tmax (e.g. MEG, where
191 // TminPsat = Tmax = 373.15 K leaves no valid range): #2209. Throw a
192 // controlled exception instead so callers see the limitation rather
193 // than a misleading zero.
194 if (T <= this->TminPsat) {
195 throw ValueError(format("Saturation pressure is not available below TminPsat=%g K (T=%g K)", this->TminPsat, T));
196 }
197 switch (p_sat.type) {
198 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
199 return poly.evaluate(p_sat.coeffs, T, x, 0, 0, Tbase, xbase);
200 case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
201 return baseExponential(p_sat, T, 0.0);
202 case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
203 return baseLogexponential(p_sat, T, 0.0);
204 case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
205 return exp(poly.evaluate(p_sat.coeffs, T, x, 0, 0, Tbase, xbase));
206 case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
207 return basePolyOffset(p_sat, T, x);
208 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
209 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?", __FILE__,
210 __LINE__, p_sat.type));
211 default:
212 throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.", __FILE__, __LINE__, p_sat.type));
213 }
214}
216double IncompressibleFluid::Tfreeze(double p, double x) {
217 switch (T_freeze.type) {
218 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
219 return poly.evaluate(T_freeze.coeffs, p, x, 0, 0, 0.0, xbase);
220 case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
221 return baseExponential(T_freeze, x, 0.0);
222 case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
223 return baseLogexponential(T_freeze, x, 0.0);
224 case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
225 return exp(poly.evaluate(T_freeze.coeffs, p, x, 0, 0, 0.0, xbase));
226 case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
227 return basePolyOffset(T_freeze, p, x);
228 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
229 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?", __FILE__,
230 __LINE__, T_freeze.type));
231 default:
232 throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.", __FILE__, __LINE__, T_freeze.type));
233 }
234}
235
236/* Below are direct calculations of the derivatives. Nothing
237 * special is going on, we simply use the polynomial class to
238 * derive the different functions with respect to temperature.
239 */
241double IncompressibleFluid::drhodTatPx(double T, double p, double x) {
242 switch (density.type) {
243 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
244 return poly.derivative(density.coeffs, T, x, 0, 0, 0, Tbase, xbase);
245 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
246 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?", __FILE__,
247 __LINE__, density.type));
248 default:
249 throw ValueError(
250 format("%s (%d): There is no predefined way to use this function type \"[%d]\" for density.", __FILE__, __LINE__, density.type));
251 }
252}
254// with respect to temperature at constant pressure and composition
255// integrated in temperature
256double IncompressibleFluid::dsdTatPxdT(double T, double p, double x) {
257 switch (specific_heat.type) {
258 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
259 return poly.integral(specific_heat.coeffs, T, x, 0, -1, 0, Tbase, xbase);
260 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
261 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?", __FILE__,
262 __LINE__, specific_heat.type));
263 default:
264 throw ValueError(
265 format("%s (%d): There is no predefined way to use this function type \"[%d]\" for entropy.", __FILE__, __LINE__, specific_heat.type));
266 }
267}
269// with respect to temperature at constant pressure and composition
270// integrated in temperature
271double IncompressibleFluid::dhdTatPxdT(double T, double p, double x) {
272 switch (specific_heat.type) {
273 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
274 return poly.integral(specific_heat.coeffs, T, x, 0, 0, 0, Tbase, xbase);
275 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
276 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?", __FILE__,
277 __LINE__, specific_heat.type));
278 default:
279 throw ValueError(
280 format("%s (%d): There is no predefined way to use this function type \"[%d]\" for entropy.", __FILE__, __LINE__, specific_heat.type));
281 }
282}
283
285
287double IncompressibleFluid::inputFromMass(double T, double x) {
288 if (this->xid == IFRAC_PURE) {
289 return _HUGE;
290 } else if (this->xid == IFRAC_MASS) {
291 return x;
292 } else {
293 throw NotImplementedError("Mass composition conversion has not been implemented.");
294 //switch (mass2input.type) {
295 // case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
296 // return poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase are defined in the correct way
297 // break;
298 // case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
299 // return baseExponential(mass2input, x, 0.0);
300 // break;
301 // case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
302 // return baseLogexponential(mass2input, x, 0.0);
303 // break;
304 // case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
305 // return exp(poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase are defined in the correct way
306 // break;
307 // case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
308 // return basePolyOffset(mass2input, T, x);
309 // break;
310 // case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
311 // throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,mass2input.type));
312 // break;
313 // default:
314 // throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.",__FILE__,__LINE__,mass2input.type));
315 // break;
316 //}
317 //return _HUGE;
318 }
319}
320
322
324double IncompressibleFluid::inputFromVolume(double T, double x) {
325 if (this->xid == IFRAC_PURE) {
326 return _HUGE;
327 } else if (this->xid == IFRAC_VOLUME) {
328 return x;
329 } else {
330 throw NotImplementedError("Volume composition conversion has not been implemented.");
331 //switch (volume2input.type) {
332 // case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
333 // return poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase are defined in the correct way
334 // break;
335 // case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
336 // return baseExponential(volume2input, x, 0.0);
337 // break;
338 // case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
339 // return baseLogexponential(volume2input, x, 0.0);
340 // break;
341 // case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
342 // return exp(poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase are defined in the correct way
343 // break;
344 // case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
345 // return basePolyOffset(volume2input, T, x);
346 // break;
347 // case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
348 // throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,volume2input.type));
349 // break;
350 // default:
351 // throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.",__FILE__,__LINE__,volume2input.type));
352 // break;
353 //}
354 //return _HUGE;
355 }
356}
357
359
361double IncompressibleFluid::inputFromMole(double T, double x) {
362 if (this->xid == IFRAC_PURE) {
363 return _HUGE;
364 } else if (this->xid == IFRAC_MOLE) {
365 return x;
366 } else {
367 throw NotImplementedError("Mole composition conversion has not been implemented.");
368 /*
369 switch (mole2input.type) {
370 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
371 return poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase are defined in the correct way
372 break;
373 case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
374 return baseExponential(mole2input, x, 0.0);
375 break;
376 case IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL:
377 return baseLogexponential(mole2input, x, 0.0);
378 break;
379 case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
380 return exp(poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase are defined in the correct way
381 break;
382 case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
383 return basePolyOffset(mole2input, T, x);
384 break;
385 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
386 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,mole2input.type));
387 break;
388 default:
389 throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.",__FILE__,__LINE__,mole2input.type));
390 break;
391 }
392 return _HUGE;
393 */
394 }
395}
396
397/* Some functions can be inverted directly, those are listed
398 * here. It is also possible to solve for other quantities, but
399 * that involves some more sophisticated processing and is not
400 * done here, but in the backend, T(h,p) for example.
401 */
403double IncompressibleFluid::T_rho(double Dmass, double p, double x) {
404 double d_raw = Dmass; // No changes needed, no reference values...
405 switch (density.type) {
406 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
407 return poly.solve_limits(density.coeffs, x, d_raw, Tmin, Tmax, 0, 0, 0, Tbase, xbase);
408 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
409 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?", __FILE__,
410 __LINE__, specific_heat.type));
411 default:
412 throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse density.", __FILE__, __LINE__,
413 specific_heat.type));
414 }
415}
417double IncompressibleFluid::T_c(double Cmass, double p, double x) {
418 double c_raw = Cmass; // No changes needed, no reference values...
419 switch (specific_heat.type) {
420 case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
421 return poly.solve_limits(specific_heat.coeffs, x, c_raw, Tmin, Tmax, 0, 0, 0, Tbase, xbase);
422 case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
423 throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?", __FILE__,
424 __LINE__, specific_heat.type));
425 default:
426 throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse specific heat.", __FILE__,
427 __LINE__, specific_heat.type));
428 }
429}
430
431/*
432 * Some more functions to provide a single implementation
433 * of important routines.
434 * We start with the check functions that can validate input
435 * in terms of pressure p, temperature T and composition x.
436 */
438
441bool IncompressibleFluid::checkT(double T, double p, double x) {
442 if (Tmin <= 0.) throw ValueError("Please specify the minimum temperature.");
443 if (Tmax <= 0.) throw ValueError("Please specify the maximum temperature.");
444 if ((Tmin > T) || (T > Tmax)) throw ValueError(format("Your temperature %f is not between %f and %f.", T, Tmin, Tmax));
445 double TF = 0.0;
446 if (T_freeze.type != IncompressibleData::INCOMPRESSIBLE_NOT_SET) TF = Tfreeze(p, x);
447 if (T < TF) throw ValueError(format("Your temperature %f is below the freezing point of %f.", T, TF));
448 return true;
449}
450
452
458bool IncompressibleFluid::checkP(double T, double p, double x) {
459 double ps = 0.0;
460 if (p_sat.type != IncompressibleData::INCOMPRESSIBLE_NOT_SET) {
461 // psat() now throws below TminPsat (#2209). For this validation
462 // path, T below TminPsat means the saturation curve is not
463 // available — there is nothing to compare p against, so skip
464 // the p >= psat check rather than propagating the throw to
465 // callers that didn't ask for psat in the first place
466 // (e.g. PropsSI("D","P",101325,"T",300,"INCOMP::DowQ") where
467 // DowQ has TminPsat > 300 K).
468 try {
469 ps = psat(T, x);
470 } catch (const ValueError&) {
471 ps = 0.0;
472 }
473 }
474 if (p < 0.0) throw ValueError(format("You cannot use negative pressures: %f < %f. ", p, 0.0));
475 if (ps > 0.0 && p < ps) throw ValueError(format("Equations are valid for liquid phase only: %f < %f (psat). ", p, ps));
476 return true;
477}
478
480
483bool IncompressibleFluid::checkX(double x) {
484 if (xmin < 0.0 || xmin > 1.0) throw ValueError("Please specify the minimum concentration between 0 and 1.");
485 if (xmax < 0.0 || xmax > 1.0) throw ValueError("Please specify the maximum concentration between 0 and 1.");
486 if ((x < xmin * (1 - INCOMP_EPSILON)) || (x > xmax * (1 + INCOMP_EPSILON))) {
487 throw ValueError(format("Your composition %g is not between %g and %g.", x, xmin, xmax));
488 }
489 return true;
490}
491
492} /* namespace CoolProp */
493
494// Testing still needs to be enhanced.
495/* Below, I try to carry out some basic tests for both 2D and 1D
496 * polynomials as well as the exponential functions for vapour
497 * pressure etc.
498 */
499#ifdef ENABLE_CATCH
500# include <math.h>
501# include <iostream>
502# include <catch2/catch_all.hpp>
503# include "Tests/TestObjects.h"
504
505Eigen::MatrixXd makeMatrix(const std::vector<double>& coefficients) {
506 //IncompressibleClass::checkCoefficients(coefficients,18);
507 std::vector<std::vector<double>> matrix;
508 std::vector<double> tmpVector;
509
510 tmpVector.clear();
511 tmpVector.push_back(coefficients[0]);
512 tmpVector.push_back(coefficients[6]);
513 tmpVector.push_back(coefficients[11]);
514 tmpVector.push_back(coefficients[15]);
515 matrix.push_back(tmpVector);
516
517 tmpVector.clear();
518 tmpVector.push_back(coefficients[1] * 100.0);
519 tmpVector.push_back(coefficients[7] * 100.0);
520 tmpVector.push_back(coefficients[12] * 100.0);
521 tmpVector.push_back(coefficients[16] * 100.0);
522 matrix.push_back(tmpVector);
523
524 tmpVector.clear();
525 tmpVector.push_back(coefficients[2] * 100.0 * 100.0);
526 tmpVector.push_back(coefficients[8] * 100.0 * 100.0);
527 tmpVector.push_back(coefficients[13] * 100.0 * 100.0);
528 tmpVector.push_back(coefficients[17] * 100.0 * 100.0);
529 matrix.push_back(tmpVector);
530
531 tmpVector.clear();
532 tmpVector.push_back(coefficients[3] * 100.0 * 100.0 * 100.0);
533 tmpVector.push_back(coefficients[9] * 100.0 * 100.0 * 100.0);
534 tmpVector.push_back(coefficients[14] * 100.0 * 100.0 * 100.0);
535 tmpVector.push_back(0.0);
536 matrix.push_back(tmpVector);
537
538 tmpVector.clear();
539 tmpVector.push_back(coefficients[4] * 100.0 * 100.0 * 100.0 * 100.0);
540 tmpVector.push_back(coefficients[10] * 100.0 * 100.0 * 100.0 * 100.0);
541 tmpVector.push_back(0.0);
542 tmpVector.push_back(0.0);
543 matrix.push_back(tmpVector);
544
545 tmpVector.clear();
546 tmpVector.push_back(coefficients[5] * 100.0 * 100.0 * 100.0 * 100.0 * 100.0);
547 tmpVector.push_back(0.0);
548 tmpVector.push_back(0.0);
549 tmpVector.push_back(0.0);
550 matrix.push_back(tmpVector);
551
552 tmpVector.clear();
553 return CoolProp::vec_to_eigen(matrix).transpose();
554}
555
556TEST_CASE("Internal consistency checks and example use cases for the incompressible fluids", "[IncompressibleFluids]") {
557 SECTION("Test case for \"SylthermXLT\" by Dow Chemicals") {
558
559 std::vector<double> cRho;
560 cRho.push_back(+1.1563685145E+03);
561 cRho.push_back(-1.0269048032E+00);
562 cRho.push_back(-9.3506079577E-07);
563 cRho.push_back(+1.0368116627E-09);
564 CoolProp::IncompressibleData density;
565 density.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
566 density.coeffs = CoolProp::vec_to_eigen(cRho);
567
568 std::vector<double> cHeat;
569 cHeat.push_back(+1.1562261074E+03);
570 cHeat.push_back(+2.0994549103E+00);
571 cHeat.push_back(+7.7175381057E-07);
572 cHeat.push_back(-3.7008444051E-20);
573 CoolProp::IncompressibleData specific_heat;
574 specific_heat.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
575 specific_heat.coeffs = CoolProp::vec_to_eigen(cHeat);
576
577 std::vector<double> cCond;
578 cCond.push_back(+1.6121957379E-01);
579 cCond.push_back(-1.3023781944E-04);
580 cCond.push_back(-1.4395238766E-07);
581 CoolProp::IncompressibleData conductivity;
582 conductivity.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
583 conductivity.coeffs = CoolProp::vec_to_eigen(cCond);
584
585 std::vector<double> cVisc;
586 cVisc.push_back(+1.0337654989E+03);
587 cVisc.push_back(-4.3322764383E+01);
588 cVisc.push_back(+1.0715062356E+01);
589 CoolProp::IncompressibleData viscosity;
590 viscosity.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL;
591 viscosity.coeffs = CoolProp::vec_to_eigen(cVisc);
592
593 CoolProp::IncompressibleFluid XLT;
594 XLT.setName("XLT");
595 XLT.setDescription("SylthermXLT");
596 XLT.setReference("Dow Chemicals data sheet");
597 XLT.setTmax(533.15);
598 XLT.setTmin(173.15);
599 XLT.setxmax(0.0);
600 XLT.setxmin(0.0);
601 XLT.setTminPsat(533.15);
602
603 XLT.setTbase(0.0);
604 XLT.setxbase(0.0);
605
607 XLT.setDensity(density);
608 XLT.setSpecificHeat(specific_heat);
609 XLT.setViscosity(viscosity);
610 XLT.setConductivity(conductivity);
611 //XLT.setPsat(parse_coefficients(fluid_json, "saturation_pressure", false));
612 //XLT.setTfreeze(parse_coefficients(fluid_json, "T_freeze", false));
613 //XLT.setVolToMass(parse_coefficients(fluid_json, "volume2mass", false));
614 //XLT.setMassToMole(parse_coefficients(fluid_json, "mass2mole", false));
615
617 //XLT.validate();
618 double acc = 0.0001;
619 double val = 0;
620 double res = 0;
621
622 // Prepare the results and compare them to the calculated values
623 double T = 273.15 + 50;
624 double p = 10e5;
625 double x = 0.0;
626
627 // Compare density
628 val = 824.4615702148608;
629 res = XLT.rho(T, p, x);
630 {
631 CAPTURE(T);
632 CAPTURE(val);
633 CAPTURE(res);
634 CHECK(check_abs(val, res, acc));
635 }
636
637 // Compare cp
638 val = 1834.7455527670554;
639 res = XLT.c(T, p, x);
640 {
641 CAPTURE(T);
642 CAPTURE(val);
643 CAPTURE(res);
644 CHECK(check_abs(val, res, acc));
645 }
646
647 // Check property functions
648 CHECK_THROWS(XLT.s(T, p, x));
649 CHECK_THROWS(XLT.h(T, p, x));
650 CHECK_THROWS(XLT.u(T, p, x));
651
652 // Compare v
653 val = 0.0008931435169681835;
654 res = XLT.visc(T, p, x);
655 {
656 CAPTURE(T);
657 CAPTURE(val);
658 CAPTURE(res);
659 CHECK(check_abs(val, res, acc));
660 }
661
662 // Compare l
663 val = 0.10410086156049088;
664 res = XLT.cond(T, p, x);
665 {
666 CAPTURE(T);
667 CAPTURE(val);
668 CAPTURE(res);
669 CHECK(check_abs(val, res, acc));
670 }
671 }
672
673 SECTION("Test case for Methanol from SecCool") {
674
675 CoolProp::IncompressibleFluid CH3OH = CoolProp::get_incompressible_fluid("MMA2");
676
677 // Prepare the results and compare them to the calculated values
678 double acc = 0.0001;
679 double T = 273.15 + 10;
680 double p = 10e5;
681 double x = 0.25;
682 double expected = 0;
683 double actual = 0;
684
685 // Compare density
686 expected = 963.2886528091547;
687 actual = CH3OH.rho(T, p, x);
688 {
689 CAPTURE(T);
690 CAPTURE(p);
691 CAPTURE(x);
692 CAPTURE(expected);
693 CAPTURE(actual);
694 CHECK(check_abs(expected, actual, acc));
695 }
696
697 // Compare cp
698 expected = 3993.9748117022423;
699 actual = CH3OH.c(T, p, x);
700 {
701 CAPTURE(T);
702 CAPTURE(p);
703 CAPTURE(x);
704 CAPTURE(expected);
705 CAPTURE(actual);
706 CHECK(check_abs(expected, actual, acc));
707 }
708
709 // Check property functions
710 CHECK_THROWS(CH3OH.s(T, p, x));
711 CHECK_THROWS(CH3OH.h(T, p, x));
712 CHECK_THROWS(CH3OH.u(T, p, x));
713
714 // Compare viscosity
715 expected = 0.0023970245009602097; // Pa-s // Melinder 3.48 mPas @ 24.87 wt%, 0C
716 actual = CH3OH.visc(T, p, x);
717 {
718 CAPTURE(T);
719 CAPTURE(p);
720 CAPTURE(x);
721 CAPTURE(expected);
722 CAPTURE(actual);
723 std::string errmsg = CoolProp::get_global_param_string("errstring");
724 CAPTURE(errmsg);
725 CHECK(check_abs(expected, actual, acc));
726 }
727
728 // Compare conductivity
729 expected = 0.44791148414693727;
730 actual = CH3OH.cond(T, p, x);
731 {
732 CAPTURE(T);
733 CAPTURE(p);
734 CAPTURE(x);
735 CAPTURE(expected);
736 CAPTURE(actual);
737 std::string errmsg = CoolProp::get_global_param_string("errstring");
738 CAPTURE(errmsg);
739 CHECK(check_abs(expected, actual, acc));
740 }
741
742 // Compare Tfreeze
743 expected = -20.02 + 273.15; // 253.1293105454671;
744 actual = CH3OH.Tfreeze(p, x);
745 {
746 CAPTURE(T);
747 CAPTURE(p);
748 CAPTURE(x);
749 CAPTURE(expected);
750 CAPTURE(actual);
751 std::string errmsg = CoolProp::get_global_param_string("errstring");
752 CAPTURE(errmsg);
753 CHECK(check_abs(expected, actual, acc));
754 }
755 }
756}
757
758#endif /* ENABLE_CATCH */
Generated by doxygen 1.9.4