CoolProp  6.6.0
An open-source fluid property and humid air property database
Ancillaries.cpp
Go to the documentation of this file.
1 #include "Ancillaries.h"
2 #include "DataStructures.h"
3 #include "AbstractState.h"
4 
5 #if defined(ENABLE_CATCH)
6 
7 # include "crossplatform_shared_ptr.h"
8 # include <catch2/catch_all.hpp>
9 
10 #endif
11 
12 namespace CoolProp {
13 
14 SaturationAncillaryFunction::SaturationAncillaryFunction(rapidjson::Value& json_code) {
15  std::string type = cpjson::get_string(json_code, "type");
16  if (!type.compare("rational_polynomial")) {
17  this->type = TYPE_RATIONAL_POLYNOMIAL;
18  num_coeffs = vec_to_eigen(cpjson::get_double_array(json_code["A"]));
19  den_coeffs = vec_to_eigen(cpjson::get_double_array(json_code["B"]));
20  max_abs_error = cpjson::get_double(json_code, "max_abs_error");
21  try {
22  Tmin = cpjson::get_double(json_code, "Tmin");
23  Tmax = cpjson::get_double(json_code, "Tmax");
24  } catch (...) {
25  Tmin = _HUGE;
26  Tmax = _HUGE;
27  }
28  } else {
29  if (!type.compare("rhoLnoexp"))
30  this->type = TYPE_NOT_EXPONENTIAL;
31  else
32  this->type = TYPE_EXPONENTIAL;
33  n = cpjson::get_double_array(json_code["n"]);
34  N = n.size();
35  s = n;
36  t = cpjson::get_double_array(json_code["t"]);
37  Tmin = cpjson::get_double(json_code, "Tmin");
38  Tmax = cpjson::get_double(json_code, "Tmax");
39  reducing_value = cpjson::get_double(json_code, "reducing_value");
40  using_tau_r = cpjson::get_bool(json_code, "using_tau_r");
41  T_r = cpjson::get_double(json_code, "T_r");
42  }
43 };
44 
45 double SaturationAncillaryFunction::evaluate(double T) {
46  if (type == TYPE_NOT_SET) {
47  throw ValueError(format("type not set"));
48  } else if (type == TYPE_RATIONAL_POLYNOMIAL) {
49  Polynomial2D poly;
50  return poly.evaluate(num_coeffs, T) / poly.evaluate(den_coeffs, T);
51  } else {
52  double THETA = 1 - T / T_r;
53 
54  for (std::size_t i = 0; i < N; ++i) {
55  s[i] = n[i] * pow(THETA, t[i]);
56  }
57  double summer = std::accumulate(s.begin(), s.end(), 0.0);
58 
59  if (type == TYPE_NOT_EXPONENTIAL) {
60  return reducing_value * (1 + summer);
61  } else {
62  double tau_r_value;
63  if (using_tau_r)
64  tau_r_value = T_r / T;
65  else
66  tau_r_value = 1.0;
67  return reducing_value * exp(tau_r_value * summer);
68  }
69  }
70 }
71 double SaturationAncillaryFunction::invert(double value, double min_bound, double max_bound) {
72  // Invert the ancillary curve to get the temperature as a function of the output variable
73  // Define the residual to be driven to zero
74  class solver_resid : public FuncWrapper1D
75  {
76  public:
77  SaturationAncillaryFunction* anc;
78  CoolPropDbl value;
79 
80  solver_resid(SaturationAncillaryFunction* anc, CoolPropDbl value) : anc(anc), value(value) {}
81 
82  double call(double T) {
83  CoolPropDbl current_value = anc->evaluate(T);
84  return current_value - value;
85  }
86  };
87  solver_resid resid(this, value);
88  std::string errstring;
89  if (min_bound < 0) {
90  min_bound = Tmin - 0.01;
91  }
92  if (max_bound < 0) {
93  max_bound = Tmax;
94  }
95 
96  try {
97  // Safe to expand the domain a little bit to lower temperature, absolutely cannot exceed Tmax
98  // because then you get (negative number)^(double) which is undefined.
99  return Brent(resid, min_bound, max_bound, DBL_EPSILON, 1e-10, 100);
100  } catch (...) {
101  return ExtrapolatingSecant(resid,max_bound, -0.01, 1e-12, 100);
102  }
103 }
104 
105 void MeltingLineVariables::set_limits(void) {
106  if (type == MELTING_LINE_SIMON_TYPE) {
107 
108  // Fill in the min and max pressures for each part
109  for (std::size_t i = 0; i < simon.parts.size(); ++i) {
110  MeltingLinePiecewiseSimonSegment& part = simon.parts[i];
111  part.p_min = part.p_0 + part.a * (pow(part.T_min / part.T_0, part.c) - 1);
112  part.p_max = part.p_0 + part.a * (pow(part.T_max / part.T_0, part.c) - 1);
113  }
114  pmin = simon.parts.front().p_min;
115  pmax = simon.parts.back().p_max;
116  Tmin = simon.parts.front().T_min;
117  Tmax = simon.parts.back().T_max;
118  } else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE) {
119  // Fill in the min and max pressures for each part
120  for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i) {
121  MeltingLinePiecewisePolynomialInTrSegment& part = polynomial_in_Tr.parts[i];
122  part.p_min = part.evaluate(part.T_min);
123  part.p_max = part.evaluate(part.T_max);
124  }
125  Tmin = polynomial_in_Tr.parts.front().T_min;
126  pmin = polynomial_in_Tr.parts.front().p_min;
127  Tmax = polynomial_in_Tr.parts.back().T_max;
128  pmax = polynomial_in_Tr.parts.back().p_max;
129  } else if (type == MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE) {
130  // Fill in the min and max pressures for each part
131  for (std::size_t i = 0; i < polynomial_in_Theta.parts.size(); ++i) {
132  MeltingLinePiecewisePolynomialInThetaSegment& part = polynomial_in_Theta.parts[i];
133  part.p_min = part.evaluate(part.T_min);
134  part.p_max = part.evaluate(part.T_max);
135  }
136  Tmin = polynomial_in_Theta.parts.front().T_min;
137  pmin = polynomial_in_Theta.parts.front().p_min;
138  Tmax = polynomial_in_Theta.parts.back().T_max;
139  pmax = polynomial_in_Theta.parts.back().p_max;
140  } else {
141  throw ValueError("only Simon supported now");
142  }
143 }
144 
145 CoolPropDbl MeltingLineVariables::evaluate(int OF, int GIVEN, CoolPropDbl value) {
146  if (type == MELTING_LINE_NOT_SET) {
147  throw ValueError("Melting line curve not set");
148  }
149  if (OF == iP_max) {
150  return pmax;
151  } else if (OF == iP_min) {
152  return pmin;
153  } else if (OF == iT_max) {
154  return Tmax;
155  } else if (OF == iT_min) {
156  return Tmin;
157  } else if (OF == iP && GIVEN == iT) {
158  CoolPropDbl T = value;
159  if (type == MELTING_LINE_SIMON_TYPE) {
160  // Need to find the right segment
161  for (std::size_t i = 0; i < simon.parts.size(); ++i) {
162  MeltingLinePiecewiseSimonSegment& part = simon.parts[i];
163  if (is_in_closed_range(part.T_min, part.T_max, T)) {
164  return part.p_0 + part.a * (pow(T / part.T_0, part.c) - 1);
165  }
166  }
167  throw ValueError("unable to calculate melting line (p,T) for Simon curve");
168  } else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE) {
169  // Need to find the right segment
170  for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i) {
171  MeltingLinePiecewisePolynomialInTrSegment& part = polynomial_in_Tr.parts[i];
172  if (is_in_closed_range(part.T_min, part.T_max, T)) {
173  return part.evaluate(T);
174  }
175  }
176  throw ValueError("unable to calculate melting line (p,T) for polynomial_in_Tr curve");
177  } else if (type == MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE) {
178  // Need to find the right segment
179  for (std::size_t i = 0; i < polynomial_in_Theta.parts.size(); ++i) {
180  MeltingLinePiecewisePolynomialInThetaSegment& part = polynomial_in_Theta.parts[i];
181  if (is_in_closed_range(part.T_min, part.T_max, T)) {
182  return part.evaluate(T);
183  }
184  }
185  throw ValueError("unable to calculate melting line (p,T) for polynomial_in_Theta curve");
186  } else {
187  throw ValueError(format("Invalid melting line type [%d]", type));
188  }
189  } else {
190  if (type == MELTING_LINE_SIMON_TYPE) {
191  // Need to find the right segment
192  for (std::size_t i = 0; i < simon.parts.size(); ++i) {
193  MeltingLinePiecewiseSimonSegment& part = simon.parts[i];
194  // p = part.p_0 + part.a*(pow(T/part.T_0,part.c)-1);
195  CoolPropDbl T = pow((value - part.p_0) / part.a + 1, 1 / part.c) * part.T_0;
196  if (T >= part.T_0 && T <= part.T_max) {
197  return T;
198  }
199  }
200  throw ValueError(format("unable to calculate melting line T(p) for Simon curve for p=%Lg; bounds are %Lg,%Lg Pa", value, pmin, pmax));
201  } else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE) {
202  class solver_resid : public FuncWrapper1D
203  {
204  public:
205  MeltingLinePiecewisePolynomialInTrSegment* part;
206  CoolPropDbl given_p;
207  solver_resid(MeltingLinePiecewisePolynomialInTrSegment* part, CoolPropDbl p) : part(part), given_p(p){};
208  double call(double T) {
209 
210  CoolPropDbl calc_p = part->evaluate(T);
211 
212  // Difference between the two is to be driven to zero
213  return given_p - calc_p;
214  };
215  };
216 
217  // Need to find the right segment
218  for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i) {
219  MeltingLinePiecewisePolynomialInTrSegment& part = polynomial_in_Tr.parts[i];
220  if (is_in_closed_range(part.p_min, part.p_max, value)) {
221  solver_resid resid(&part, value);
222  double T = Brent(resid, part.T_min, part.T_max, DBL_EPSILON, 1e-12, 100);
223  return T;
224  }
225  }
226  throw ValueError(
227  format("unable to calculate melting line T(p) for polynomial_in_Theta curve for p=%Lg; bounds are %Lg,%Lg Pa", value, pmin, pmax));
228  } else if (type == MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE) {
229 
230  class solver_resid : public FuncWrapper1D
231  {
232  public:
233  MeltingLinePiecewisePolynomialInThetaSegment* part;
234  CoolPropDbl given_p;
235  solver_resid(MeltingLinePiecewisePolynomialInThetaSegment* part, CoolPropDbl p) : part(part), given_p(p){};
236  double call(double T) {
237 
238  CoolPropDbl calc_p = part->evaluate(T);
239 
240  // Difference between the two is to be driven to zero
241  return given_p - calc_p;
242  };
243  };
244 
245  // Need to find the right segment
246  for (std::size_t i = 0; i < polynomial_in_Theta.parts.size(); ++i) {
247  MeltingLinePiecewisePolynomialInThetaSegment& part = polynomial_in_Theta.parts[i];
248  if (is_in_closed_range(part.p_min, part.p_max, value)) {
249  solver_resid resid(&part, value);
250  double T = Brent(resid, part.T_min, part.T_max, DBL_EPSILON, 1e-12, 100);
251  return T;
252  }
253  }
254 
255  throw ValueError(
256  format("unable to calculate melting line T(p) for polynomial_in_Theta curve for p=%Lg; bounds are %Lg,%Lg Pa", value, pmin, pmax));
257  } else {
258  throw ValueError(format("Invalid melting line type T(p) [%d]", type));
259  }
260  }
261 }
262 
263 }; /* namespace CoolProp */
264 
265 #if defined(ENABLE_CATCH)
266 TEST_CASE("Water melting line", "[melting]") {
267  shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS", "water"));
268  int iT = CoolProp::iT, iP = CoolProp::iP;
269  SECTION("Ice Ih-liquid") {
270  double actual = AS->melting_line(iT, iP, 138.268e6);
271  double expected = 260.0;
272  CAPTURE(actual);
273  CAPTURE(expected);
274  CHECK(std::abs(actual - expected) < 0.01);
275  }
276  SECTION("Ice III-liquid") {
277  double actual = AS->melting_line(iT, iP, 268.685e6);
278  double expected = 254;
279  CAPTURE(actual);
280  CAPTURE(expected);
281  CHECK(std::abs(actual - expected) < 0.01);
282  }
283  SECTION("Ice V-liquid") {
284  double actual = AS->melting_line(iT, iP, 479.640e6);
285  double expected = 265;
286  CAPTURE(actual);
287  CAPTURE(expected);
288  CHECK(std::abs(actual - expected) < 0.01);
289  }
290  SECTION("Ice VI-liquid") {
291  double actual = AS->melting_line(iT, iP, 1356.76e6);
292  double expected = 320;
293  CAPTURE(actual);
294  CAPTURE(expected);
295  CHECK(std::abs(actual - expected) < 1);
296  }
297 }
298 
299 TEST_CASE("Tests for values from melting lines", "[melting]") {
300  int iT = CoolProp::iT, iP = CoolProp::iP;
301  std::vector<std::string> fluids = strsplit(CoolProp::get_global_param_string("fluids_list"), ',');
302  for (std::size_t i = 0; i < fluids.size(); ++i) {
303  shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS", fluids[i]));
304 
305  // Water has its own better tests; skip fluids without melting line
306  if (!AS->has_melting_line() || !fluids[i].compare("Water")) {
307  continue;
308  }
309 
310  double pmax = AS->melting_line(CoolProp::iP_max, iT, 0);
311  double pmin = AS->melting_line(CoolProp::iP_min, iT, 0);
312  double Tmax = AS->melting_line(CoolProp::iT_max, iT, 0);
313  double Tmin = AS->melting_line(CoolProp::iT_min, iT, 0);
314 
315  // See https://groups.google.com/forum/?fromgroups#!topic/catch-forum/mRBKqtTrITU
316  std::ostringstream ss0;
317  ss0 << "Check melting line limits for fluid " << fluids[i];
318  SECTION(ss0.str(), "") {
319  CAPTURE(Tmin);
320  CAPTURE(Tmax);
321  CAPTURE(pmin);
322  CAPTURE(pmax);
323  CHECK(Tmax > Tmin);
324  CHECK(pmax > pmin);
325  CHECK(pmin > 0);
326  }
327  for (double p = 0.1 * (pmax - pmin) + pmin; p < pmax; p += 0.2 * (pmax - pmin)) {
328  // See https://groups.google.com/forum/?fromgroups#!topic/catch-forum/mRBKqtTrITU
329  std::ostringstream ss1;
330  ss1 << "Melting line for " << fluids[i] << " at p=" << p;
331  SECTION(ss1.str(), "") {
332  double actual_T = AS->melting_line(iT, iP, p);
333  CAPTURE(Tmin);
334  CAPTURE(Tmax);
335  CAPTURE(actual_T);
336  CHECK(actual_T > Tmin);
337  CHECK(actual_T < Tmax);
338  }
339  }
340  // See https://groups.google.com/forum/?fromgroups#!topic/catch-forum/mRBKqtTrITU
341  std::ostringstream ss2;
342  ss2 << "Ensure melting line valid for " << fluids[i] << " @ EOS pmax";
343  SECTION(ss2.str(), "") {
344  double actual_T;
345  double EOS_pmax = AS->pmax();
346  CAPTURE(EOS_pmax);
347  CHECK_NOTHROW(actual_T = AS->melting_line(iT, iP, EOS_pmax));
348  CAPTURE(actual_T);
349  }
350  }
351 }
352 
353 TEST_CASE("Test that hs_anchor enthalpy/entropy agrees with EOS", "[ancillaries]") {
354  std::vector<std::string> fluids = strsplit(CoolProp::get_global_param_string("fluids_list"), ',');
355  for (std::size_t i = 0; i < fluids.size(); ++i) {
356  shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS", fluids[i]));
357 
358  CoolProp::SimpleState hs_anchor = AS->get_state("hs_anchor");
359 
360  // See https://groups.google.com/forum/?fromgroups#!topic/catch-forum/mRBKqtTrITU
361  std::ostringstream ss1;
362  ss1 << "Check hs_anchor for " << fluids[i];
363  SECTION(ss1.str(), "") {
364  std::string note =
365  "The enthalpy and entropy are hardcoded in the fluid JSON files. They MUST agree with the values calculated by the EOS";
366  AS->update(CoolProp::DmolarT_INPUTS, hs_anchor.rhomolar, hs_anchor.T);
367  double EOS_hmolar = AS->hmolar();
368  double EOS_smolar = AS->smolar();
369  CAPTURE(hs_anchor.hmolar);
370  CAPTURE(hs_anchor.smolar);
371  CAPTURE(EOS_hmolar);
372  CAPTURE(EOS_smolar);
373  CHECK(std::abs(EOS_hmolar - hs_anchor.hmolar) < 1e-3);
374  CHECK(std::abs(EOS_smolar - hs_anchor.smolar) < 1e-3);
375  }
376  }
377 }
378 
379 TEST_CASE("Surface tension", "[surface_tension]") {
380  SECTION("from PropsSI") {
381  CHECK(ValidNumber(CoolProp::PropsSI("surface_tension", "T", 300, "Q", 0, "Water")));
382  }
383  SECTION("from saturation_ancillary") {
384  CHECK(ValidNumber(CoolProp::saturation_ancillary("Water", "surface_tension", 0, "T", 300)));
385  }
386  SECTION("from AbstractState") {
387  shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS", "Water"));
388  AS->update(CoolProp::QT_INPUTS, 0, 300);
389  CHECK_NOTHROW(AS->surface_tension());
390  }
391 }
392 #endif
Generated on Wed Nov 29 2023 14:26:24 for CoolProp by doxygen 1.9.1