https://mooseframework.inl.gov
Loading...
Searching...
No Matches
SodiumSaturationFluidProperties.C
Go to the documentation of this file.
1//* This file is part of the MOOSE framework
2//* https://mooseframework.inl.gov
3//*
4//* All rights reserved, see COPYRIGHT for full restrictions
5//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6//*
7//* Licensed under LGPL 2.1, please see LICENSE for details
8//* https://www.gnu.org/licenses/lgpl-2.1.html
9
11#include "NewtonInversion.h"
12
14
17{
19 params.addClassDescription("Fluid properties for liquid sodium at saturation conditions");
20 return params;
21}
22
27
28std::string
30{
31 return "sodium_sat";
32}
33
34Real
36{
37 return 22.989769E-3;
38}
39
40Real
41SodiumSaturationFluidProperties::rho_from_p_T(Real /* pressure */, Real temperature) const
42{
43 return 1.00423e3 - 0.21390 * temperature - 1.1046e-5 * temperature * temperature;
44}
45
46void
48 Real pressure, Real temperature, Real & rho, Real & drho_dp, Real & drho_dT) const
49{
50 rho = rho_from_p_T(pressure, temperature);
51 drho_dp = 0.0;
52 drho_dT = -0.21390 - 1.1046e-5 * 2 * temperature;
53}
54
55void
57 const ADReal & temperature,
58 ADReal & rho,
59 ADReal & drho_dp,
60 ADReal & drho_dT) const
61{
62 rho = SinglePhaseFluidProperties::rho_from_p_T(pressure, temperature);
63 drho_dp = 0.0;
64 drho_dT = -0.21390 - 1.1046e-5 * 2 * temperature;
65}
66
67Real
68SodiumSaturationFluidProperties::rho_from_p_s(Real pressure, Real entropy) const
69{
70 auto entropy_from_p_T = [&](Real p, Real T, Real & s, Real & ds_dp, Real & ds_dT)
71 { s_from_p_T(p, T, s, ds_dp, ds_dT); };
72 const Real temperature = FluidPropertiesUtils::NewtonSolve(pressure,
73 entropy,
76 entropy_from_p_T,
77 name() + "::rho_from_p_s",
80 .first;
81 return rho_from_p_T(pressure, temperature);
82}
83
84void
86 Real pressure, Real entropy, Real & rho, Real & drho_dp, Real & drho_ds) const
87{
88 auto entropy_from_p_T = [&](Real p, Real T, Real & s, Real & ds_dp, Real & ds_dT)
89 { s_from_p_T(p, T, s, ds_dp, ds_dT); };
90 const Real temperature = FluidPropertiesUtils::NewtonSolve(pressure,
91 entropy,
94 entropy_from_p_T,
95 name() + "::rho_from_p_s",
98 .first;
99
100 Real s, ds_dp, ds_dT;
101 s_from_p_T(pressure, temperature, s, ds_dp, ds_dT);
102 const Real dT_dp = -ds_dp / ds_dT;
103 const Real dT_ds = 1 / ds_dT;
104
105 Real drho_dp_T, drho_dT;
106 rho_from_p_T(pressure, temperature, rho, drho_dp_T, drho_dT);
107 drho_dp = drho_dp_T + drho_dT * dT_dp;
108 drho_ds = drho_dT * dT_ds;
109}
110
111Real
112SodiumSaturationFluidProperties::v_from_p_T(Real pressure, Real temperature) const
113{
114 return 1.0 / rho_from_p_T(pressure, temperature);
115}
116
117void
119 Real temperature, Real & v, Real & dv_dT, Real & d2v_dT2, Real & d3v_dT3) const
120{
121 const Real rho = rho_from_p_T(_reference_pressure, temperature);
122 const Real drho_dT = -0.21390 - 2 * 1.1046e-5 * temperature;
123 const Real d2rho_dT2 = -2 * 1.1046e-5;
124
125 v = 1.0 / rho;
126 dv_dT = -drho_dT / (rho * rho);
127 d2v_dT2 = 2 * drho_dT * drho_dT / (rho * rho * rho) - d2rho_dT2 / (rho * rho);
128 d3v_dT3 = 6 * drho_dT * d2rho_dT2 / (rho * rho * rho) -
129 6 * drho_dT * drho_dT * drho_dT / (rho * rho * rho * rho);
130}
131
132Real
134{
135 const Real t2 = temperature * temperature;
136 return 3.7782E-10 * t2 * t2 - 1.7191E-6 * t2 * temperature + 3.0921E-3 * t2 -
137 2.4560 * temperature + 1972.0;
138}
139
140Real
142{
143 const Real t2 = temperature * temperature;
144 return 4 * 3.7782E-10 * t2 * temperature - 3 * 1.7191E-6 * t2 + 2 * 3.0921e-3 * temperature -
145 2.456;
146}
147
148void
150 Real pressure, Real temperature, Real & v, Real & dv_dp, Real & dv_dT) const
151{
152 v = v_from_p_T(pressure, temperature);
153 dv_dp = 0.0;
154
155 Real drho_dT = -0.21390 - 1.1046e-5 * 2 * temperature;
156 dv_dT = -v * v * drho_dT;
157}
158
159Real
161{
162 Real temperature = T_from_v_e(v, e);
163 Real dv_dT, d2v_dT2, d3v_dT3;
164 specific_volume_derivatives(temperature, v, dv_dT, d2v_dT2, d3v_dT3);
165
166 const Real h = h_from_p_T(_reference_pressure, temperature);
167 return _reference_pressure + (h - _reference_pressure * v - e) / (temperature * dv_dT);
168}
169
170void
172 Real v, Real e, Real & pressure, Real & dp_dv, Real & dp_de) const
173{
174 Real temperature, dT_dv, dT_de;
175 T_from_v_e(v, e, temperature, dT_dv, dT_de);
176
177 Real dv_dT, d2v_dT2, d3v_dT3;
178 specific_volume_derivatives(temperature, v, dv_dT, d2v_dT2, d3v_dT3);
179
180 const Real h = h_from_p_T(_reference_pressure, temperature);
181 const Real cp = cp0_from_T(temperature);
182 const Real numerator = h - _reference_pressure * v - e;
183 const Real denominator = temperature * dv_dT;
184 const Real dnumerator_dv = (cp - _reference_pressure * dv_dT) * dT_dv;
185 const Real ddenominator_dv = (dv_dT + temperature * d2v_dT2) * dT_dv;
186
187 pressure = _reference_pressure + numerator / denominator;
188 dp_dv = (dnumerator_dv * denominator - numerator * ddenominator_dv) / (denominator * denominator);
189 dp_de = -1.0 / denominator;
190}
191
192Real
194{
195 // From inversion of second order polynomial form of rho(T)
196 mooseAssert(0.2139 * 0.2139 + 4 * 1.1046e5 * (1.00423e3 - 1 / v) > 0,
197 "Specific volume out of bounds");
198 return (0.2139 - std::sqrt(0.2139 * 0.2139 + 4 * 1.1046e-5 * (1.00423e3 - 1 / v))) /
199 (2 * -1.1046e-5);
200}
201
202void
204 Real v, Real e, Real & temperature, Real & dT_dv, Real & dT_de) const
205{
206 temperature = T_from_v_e(v, e);
207 const Real drho_dT = -0.21390 - 1.1046e-5 * 2 * temperature;
208 dT_dv = -1 / (v * v * drho_dT);
209 dT_de = 0;
210}
211
212Real
214{
215 const Real temperature = T_from_v_e(v, e);
216 // Equation (3) in the saturated-sodium speed-of-sound section of Fink and Leibowitz (1979),
217 // valid from 370.98 K to 1173 K.
218 return 2660.7 - 0.37667 * temperature - 9.0356e-5 * temperature * temperature;
219}
220
221void
223 Real v, Real e, Real & c, Real & dc_dv, Real & dc_de) const
224{
225 Real temperature, dT_dv, dT_de;
226 T_from_v_e(v, e, temperature, dT_dv, dT_de);
227
228 c = c_from_v_e(v, e);
229 const Real dc_dT = -0.37667 - 2 * 9.0356e-5 * temperature;
230 dc_dv = dc_dT * dT_dv;
231 dc_de = dc_dT * dT_de;
232}
233
234Real
235SodiumSaturationFluidProperties::h_from_p_T(Real pressure, Real temperature) const
236{
237 Real t2 = temperature * temperature;
238 const Real h0 = 3.7782E-10 * t2 * t2 * temperature / 5 - 1.7191E-6 * t2 * t2 / 4.0 +
239 3.0921E-3 * t2 * temperature / 3.0 - 2.4560 * t2 / 2.0 + 1972.0 * temperature -
240 401088.7;
241
242 Real v, dv_dT, d2v_dT2, d3v_dT3;
243 specific_volume_derivatives(temperature, v, dv_dT, d2v_dT2, d3v_dT3);
244 return h0 + (pressure - _reference_pressure) * (v - temperature * dv_dT);
245}
246
247void
249 Real pressure, Real temperature, Real & h, Real & dh_dp, Real & dh_dT) const
250{
251 h = h_from_p_T(pressure, temperature);
252
253 Real v, dv_dT, d2v_dT2, d3v_dT3;
254 specific_volume_derivatives(temperature, v, dv_dT, d2v_dT2, d3v_dT3);
255 dh_dp = v - temperature * dv_dT;
256 dh_dT = cp0_from_T(temperature) - (pressure - _reference_pressure) * temperature * d2v_dT2;
257}
258
259Real
260SodiumSaturationFluidProperties::T_from_p_h(Real pressure, Real enthalpy) const
261{
262 auto enthalpy_from_p_T = [&](Real p, Real T, Real & h, Real & dh_dp, Real & dh_dT)
263 { h_from_p_T(p, T, h, dh_dp, dh_dT); };
264
265 return FluidPropertiesUtils::NewtonSolve(pressure,
266 enthalpy,
269 enthalpy_from_p_T,
270 name() + "::T_from_p_h",
273 .first;
274}
275
276void
278 Real pressure, Real enthalpy, Real & temperature, Real & dT_dp, Real & dT_dh) const
279{
280 temperature = T_from_p_h(pressure, enthalpy);
281
282 Real h, dh_dp, dh_dT;
283 h_from_p_T(pressure, temperature, h, dh_dp, dh_dT);
284 dT_dp = -dh_dp / dh_dT;
285 dT_dh = 1.0 / dh_dT;
286}
287
288Real
289SodiumSaturationFluidProperties::s_from_p_T(Real pressure, Real temperature) const
290{
291 // Use a first-order pressure extension consistent with the Gibbs relation for v=v(T).
292 constexpr Real reference_temperature = 370.98;
293 const auto temperature_part = [](Real T)
294 {
295 const Real T2 = T * T;
296 return 3.7782e-10 * T2 * T2 / 4.0 - 1.7191e-6 * T2 * T / 3.0 + 3.0921e-3 * T2 / 2.0 -
297 2.4560 * T + 1972.0 * std::log(T);
298 };
299
300 Real v, dv_dT, d2v_dT2, d3v_dT3;
301 specific_volume_derivatives(temperature, v, dv_dT, d2v_dT2, d3v_dT3);
302 return temperature_part(temperature) - temperature_part(reference_temperature) -
303 (pressure - _reference_pressure) * dv_dT;
304}
305
306void
308 Real pressure, Real temperature, Real & s, Real & ds_dp, Real & ds_dT) const
309{
310 s = s_from_p_T(pressure, temperature);
311
312 Real v, dv_dT, d2v_dT2, d3v_dT3;
313 specific_volume_derivatives(temperature, v, dv_dT, d2v_dT2, d3v_dT3);
314 ds_dp = -dv_dT;
315 ds_dT = cp_from_p_T(pressure, temperature) / temperature;
316}
317
318Real
320{
321 return s_from_p_T(p_from_v_e(v, e), T_from_v_e(v, e));
322}
323
324void
326 Real v, Real e, Real & s, Real & ds_dv, Real & ds_de) const
327{
328 Real pressure, dp_dv, dp_de;
329 p_from_v_e(v, e, pressure, dp_dv, dp_de);
330
331 Real temperature, dT_dv, dT_de;
332 T_from_v_e(v, e, temperature, dT_dv, dT_de);
333
334 Real ds_dp, ds_dT;
335 s_from_p_T(pressure, temperature, s, ds_dp, ds_dT);
336 ds_dv = ds_dp * dp_dv + ds_dT * dT_dv;
337 ds_de = ds_dp * dp_de + ds_dT * dT_de;
338}
339
340Real
341SodiumSaturationFluidProperties::e_from_p_T(Real pressure, Real temperature) const
342{
343 // definition of h = e + p * v
344 Real v = v_from_p_T(pressure, temperature);
345 Real h = h_from_p_T(pressure, temperature);
346 return h - pressure * v;
347}
348
349void
351 Real pressure, Real temperature, Real & e, Real & de_dp, Real & de_dT) const
352{
353 Real v, dv_dp, dv_dT;
354 v_from_p_T(pressure, temperature, v, dv_dp, dv_dT);
355 Real h, dh_dp, dh_dT;
356 h_from_p_T(pressure, temperature, h, dh_dp, dh_dT);
357 e = h - pressure * v;
358
359 // definition of e = h - p * v
360 de_dp = dh_dp - v - pressure * dv_dp;
361
362 de_dT = dh_dT - pressure * dv_dT;
363}
364
365Real
367{
368 return e_from_p_T(pressure, T_from_v_e(1 / rho, 0));
369}
370
371void
373 Real pressure, Real rho, Real & e, Real & de_dp, Real & de_drho) const
374{
375 const Real v = 1 / rho;
376 Real temperature, dT_dv, dT_de;
377 T_from_v_e(v, 0, temperature, dT_dv, dT_de);
378
379 Real de_dp_T, de_dT;
380 e_from_p_T(pressure, temperature, e, de_dp_T, de_dT);
381 de_dp = de_dp_T;
382 de_drho = de_dT * dT_dv * -v * v;
383}
384
385Real
386SodiumSaturationFluidProperties::cp_from_p_T(Real pressure, Real temperature) const
387{
388 Real v, dv_dT, d2v_dT2, d3v_dT3;
389 specific_volume_derivatives(temperature, v, dv_dT, d2v_dT2, d3v_dT3);
390 return cp0_from_T(temperature) - (pressure - _reference_pressure) * temperature * d2v_dT2;
391}
392
393void
395 Real pressure, Real temperature, Real & cp, Real & dcp_dp, Real & dcp_dT) const
396{
397 cp = cp_from_p_T(pressure, temperature);
398 Real v, dv_dT, d2v_dT2, d3v_dT3;
399 specific_volume_derivatives(temperature, v, dv_dT, d2v_dT2, d3v_dT3);
400 dcp_dp = -temperature * d2v_dT2;
401 dcp_dT = dcp0_dT_from_T(temperature) -
402 (pressure - _reference_pressure) * (d2v_dT2 + temperature * d3v_dT3);
403}
404
405Real
407{
408 return cp_from_p_T(p_from_v_e(v, e), T_from_v_e(v, e));
409}
410
411void
413 Real v, Real e, Real & cp, Real & dcp_dv, Real & dcp_de) const
414{
415 Real pressure, dp_dv, dp_de;
416 p_from_v_e(v, e, pressure, dp_dv, dp_de);
417 Real temperature, dT_dv, dT_de;
418 T_from_v_e(v, e, temperature, dT_dv, dT_de);
419
420 Real dcp_dp, dcp_dT;
421 cp_from_p_T(pressure, temperature, cp, dcp_dp, dcp_dT);
422 dcp_dv = dcp_dp * dp_dv + dcp_dT * dT_dv;
423 dcp_de = dcp_dp * dp_de + dcp_dT * dT_de;
424}
425
426Real
427SodiumSaturationFluidProperties::cv_from_p_T(Real /* pressure */, Real temperature) const
428{
429 Real t2 = temperature * temperature;
430 return 1.0369E-8 * temperature * t2 + 3.7164E-4 * t2 - 1.0494 * temperature + 1582.6;
431}
432
433void
435 Real pressure, Real temperature, Real & cv, Real & dcv_dp, Real & dcv_dT) const
436{
437 cv = cv_from_p_T(pressure, temperature);
438 dcv_dp = 0.0;
439 dcv_dT = 3 * 1.0369e-8 * temperature * temperature + 2 * 3.7164e-4 * temperature - 1.0494;
440}
441
442Real
444{
445 return cv_from_p_T(p_from_v_e(v, e), T_from_v_e(v, e));
446}
447
448void
450 Real v, Real e, Real & cv, Real & dcv_dv, Real & dcv_de) const
451{
452 Real pressure, dp_dv, dp_de;
453 p_from_v_e(v, e, pressure, dp_dv, dp_de);
454 Real temperature, dT_dv, dT_de;
455 T_from_v_e(v, e, temperature, dT_dv, dT_de);
456
457 Real dcv_dp, dcv_dT;
458 cv_from_p_T(pressure, temperature, cv, dcv_dp, dcv_dT);
459 dcv_dv = dcv_dp * dp_dv + dcv_dT * dT_dv;
460 dcv_de = dcv_dp * dp_de + dcv_dT * dT_de;
461}
462
463Real
464SodiumSaturationFluidProperties::mu_from_p_T(Real /*pressure*/, Real temperature) const
465{
466 return 3.6522E-5 + 0.16626 / temperature - 4.56877e1 / (temperature * temperature) +
467 2.8733E4 / (temperature * temperature * temperature);
468}
469
470void
472 Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const
473{
474 mu = this->mu_from_p_T(pressure, temperature);
475 dmu_dp = 0.0;
476
477 Real t2 = temperature * temperature;
478 dmu_dT = 0.16626 * -1 / t2 - 4.56877E1 * -2 / (temperature * t2) + 2.8733E4 * -3 / (t2 * t2);
479}
480
481Real
483{
484 return mu_from_p_T(p_from_v_e(v, e), T_from_v_e(v, e));
485}
486
487void
489 Real v, Real e, Real & mu, Real & dmu_dv, Real & dmu_de) const
490{
491 Real pressure, dp_dv, dp_de;
492 p_from_v_e(v, e, pressure, dp_dv, dp_de);
493 Real temperature, dT_dv, dT_de;
494 T_from_v_e(v, e, temperature, dT_dv, dT_de);
495
496 Real dmu_dp, dmu_dT;
497 mu_from_p_T(pressure, temperature, mu, dmu_dp, dmu_dT);
498 dmu_dv = dmu_dp * dp_dv + dmu_dT * dT_dv;
499 dmu_de = dmu_dp * dp_de + dmu_dT * dT_de;
500}
501
502Real
503SodiumSaturationFluidProperties::k_from_p_T(Real /*pressure*/, Real temperature) const
504{
505 return 1.1045e2 - 6.5112e-2 * temperature + 1.5430e-5 * temperature * temperature -
506 2.4617e-9 * temperature * temperature * temperature;
507}
508
509void
511 Real pressure, Real temperature, Real & k, Real & dk_dp, Real & dk_dT) const
512{
513 k = this->k_from_p_T(pressure, temperature);
514 dk_dp = 0.0;
515 dk_dT = -6.5112e-2 + 2 * 1.5430e-5 * temperature - 3 * 2.4617e-9 * temperature * temperature;
516}
517
518Real
520{
521 return k_from_p_T(p_from_v_e(v, e), T_from_v_e(v, e));
522}
523
524void
526 Real v, Real e, Real & k, Real & dk_dv, Real & dk_de) const
527{
528 Real pressure, dp_dv, dp_de;
529 p_from_v_e(v, e, pressure, dp_dv, dp_de);
530 Real temperature, dT_dv, dT_de;
531 T_from_v_e(v, e, temperature, dT_dv, dT_de);
532
533 Real dk_dp, dk_dT;
534 k_from_p_T(pressure, temperature, k, dk_dp, dk_dT);
535 dk_dv = dk_dp * dp_dv + dk_dT * dT_dv;
536 dk_de = dk_dp * dp_de + dk_dT * dT_de;
537}
DualNumber< Real, DNDerivativeType, true > ADReal
const double mu
const double rho
const double v
registerMooseObject("FluidPropertiesApp", SodiumSaturationFluidProperties)
void addClassDescription(const std::string &doc_string)
const std::string & name() const
Common class for single phase fluid properties.
const Real _T_initial_guess
Initial guess for temperature (or temperature used to compute the initial guess)
const bool _verbose_newton
Whether to output information about newton solves to console.
static InputParameters validParams()
const Real _tolerance
Newton's method may be used to convert between variable sets.
const unsigned int _max_newton_its
Maximum number of iterations for the variable conversion newton solves.
e e e e s T T T T T rho v v T e h
Fluid properties for liquid sodium at saturation conditions {fink}.
virtual Real v_from_p_T(Real p, Real T) const override
SodiumSaturationFluidProperties(const InputParameters &parameters)
virtual Real mu_from_p_T(Real p, Real T) const override
virtual Real e_from_p_rho(Real p, Real rho) const override
virtual Real k_from_v_e(Real v, Real e) const override
virtual Real k_from_p_T(Real p, Real T) const override
void specific_volume_derivatives(Real temperature, Real &v, Real &dv_dT, Real &d2v_dT2, Real &d3v_dT3) const
virtual Real T_from_v_e(Real v, Real e) const override
virtual Real T_from_p_h(Real p, Real h) const override
virtual Real cv_from_p_T(Real p, Real T) const override
virtual Real molarMass() const override
Molar mass [kg/mol].
virtual Real rho_from_p_T(Real p, Real T) const override
virtual Real mu_from_v_e(Real v, Real e) const override
virtual Real p_from_v_e(Real v, Real e) const override
virtual Real cp_from_v_e(Real v, Real e) const override
virtual Real cp_from_p_T(Real p, Real T) const override
virtual Real s_from_p_T(Real p, Real T) const override
virtual Real h_from_p_T(Real p, Real T) const override
virtual Real c_from_v_e(Real v, Real e) const override
virtual Real cv_from_v_e(Real v, Real e) const override
virtual Real e_from_p_T(Real p, Real T) const override
virtual std::string fluidName() const override
Fluid name.
virtual Real rho_from_p_s(Real p, Real s) const override
virtual Real s_from_v_e(Real v, Real e) const override
std::pair< T, T > NewtonSolve(const T &x, const T &y, const Real z_initial_guess, const Real tolerance, const Functor &y_from_x_z, const std::string &caller_name, const unsigned int max_its=100, const bool verbose=false)
NewtonSolve does a 1D Newton Solve to solve the equation y = f(x, z) for variable z.