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