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