CoolProp 6.8.1dev
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
12namespace CoolProp {
13
14SaturationAncillaryFunction::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
45double 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}
71double 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
105void 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
145CoolPropDbl 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)
266TEST_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
299TEST_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 = -_HUGE;
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
353TEST_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
379TEST_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 Jun 4 2025 04:12:26 for CoolProp by doxygen 1.9.4