https://mooseframework.inl.gov
IdealRealGasMixtureFluidProperties.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 
12 #include "BrentsMethod.h"
13 #include <numeric>
14 
16 
21 #define define_from_v_e_using_T_v(prop) \
22  Real IdealRealGasMixtureFluidProperties::prop##_from_v_e( \
23  Real v, Real e, const std::vector<Real> & x) const \
24  { \
25  const Real T = T_from_v_e(v, e, x); \
26  return prop##_from_T_v(T, v, x); \
27  } \
28  \
29  void IdealRealGasMixtureFluidProperties::prop##_from_v_e(Real v, \
30  Real e, \
31  const std::vector<Real> & x, \
32  Real & y, \
33  Real & dy_dv, \
34  Real & dy_de, \
35  std::vector<Real> & dy_dx) const \
36  { \
37  Real T, dT_dv_e, dT_de_v; \
38  std::vector<Real> dT_dx_ve(_n_secondary_vapors); \
39  T_from_v_e(v, e, x, T, dT_dv_e, dT_de_v, dT_dx_ve); \
40  \
41  Real dy_dT_v, dy_dv_T; \
42  std::vector<Real> dy_dx_Tv; \
43  prop##_from_T_v(T, v, x, y, dy_dT_v, dy_dv_T, dy_dx_Tv); \
44  \
45  dy_dv = dy_dv_T + dy_dT_v * dT_dv_e; \
46  dy_de = dy_dT_v * dT_de_v; \
47  dy_dx.resize(_n_secondary_vapors); \
48  for (unsigned int i = 0; i < _n_secondary_vapors; i++) \
49  dy_dx[i] = dy_dx_Tv[i] + dy_dT_v * dT_dx_ve[i]; \
50  }
51 
56 #define define_from_p_T_using_T_v(prop) \
57  Real IdealRealGasMixtureFluidProperties::prop##_from_p_T( \
58  Real p, Real T, const std::vector<Real> & x) const \
59  { \
60  const Real v = v_from_p_T(p, T, x); \
61  return prop##_from_T_v(T, v, x); \
62  } \
63  \
64  void IdealRealGasMixtureFluidProperties::prop##_from_p_T(Real p, \
65  Real T, \
66  const std::vector<Real> & x, \
67  Real & y, \
68  Real & dy_dp, \
69  Real & dy_dT, \
70  std::vector<Real> & dy_dx) const \
71  { \
72  Real v, dv_dp_T, dv_dT_p; \
73  std::vector<Real> dv_dx_pT; \
74  v_from_p_T(p, T, x, v, dv_dp_T, dv_dT_p, dv_dx_pT); \
75  \
76  Real dy_dT_v, dy_dv_T; \
77  std::vector<Real> dy_dx_Tv; \
78  prop##_from_T_v(T, v, x, y, dy_dT_v, dy_dv_T, dy_dx_Tv); \
79  \
80  dy_dp = dy_dv_T * dv_dp_T; \
81  dy_dT = dy_dT_v + dy_dv_T * dv_dT_p; \
82  dy_dx.resize(_n_secondary_vapors); \
83  for (unsigned int i = 0; i < _n_secondary_vapors; i++) \
84  dy_dx[i] = dy_dx_Tv[i] + dy_dv_T * dv_dx_pT[i]; \
85  }
86 
90 #define define_mass_specific_prop_from_T_v(prop) \
91  Real IdealRealGasMixtureFluidProperties::prop##_from_T_v( \
92  Real T, Real v, const std::vector<Real> & x) const \
93  { \
94  const Real x_primary = primaryMassFraction(x); \
95  Real y = x_primary * _fp_primary->prop##_from_T_v(T, v / x_primary); \
96  \
97  for (unsigned int i = 0; i < _n_secondary_vapors; i++) \
98  y += x[i] * _fp_secondary[i]->prop##_from_T_v(T, v / x[i]); \
99  \
100  return y; \
101  } \
102  \
103  void IdealRealGasMixtureFluidProperties::prop##_from_T_v(Real T, \
104  Real v, \
105  const std::vector<Real> & x, \
106  Real & y, \
107  Real & dy_dT, \
108  Real & dy_dv, \
109  std::vector<Real> & dy_dx) const \
110  { \
111  const Real x_primary = primaryMassFraction(x); \
112  mooseAssert(!MooseUtils::absoluteFuzzyEqual(x_primary, 0.0), "Mass fraction may not be zero"); \
113  \
114  Real y_primary, dy_dT_primary, dy_dv_primary; \
115  _fp_primary->prop##_from_T_v(T, v / x_primary, y_primary, dy_dT_primary, dy_dv_primary); \
116  y = x_primary * y_primary; \
117  dy_dT = x_primary * dy_dT_primary; \
118  dy_dv = dy_dv_primary; \
119  \
120  Real dy_dT_sec; \
121  std::vector<Real> y_sec(_n_secondary_vapors), dy_dv_sec(_n_secondary_vapors); \
122  for (unsigned int i = 0; i < _n_secondary_vapors; i++) \
123  { \
124  mooseAssert(!MooseUtils::absoluteFuzzyEqual(x[i], 0.0), "Mass fraction may not be zero"); \
125  _fp_secondary[i]->prop##_from_T_v(T, v / x[i], y_sec[i], dy_dT_sec, dy_dv_sec[i]); \
126  y += x[i] * y_sec[i]; \
127  dy_dT += x[i] * dy_dT_sec; \
128  dy_dv += dy_dv_sec[i]; \
129  } \
130  \
131  dy_dx.resize(_n_secondary_vapors); \
132  for (unsigned int i = 0; i < _n_secondary_vapors; i++) \
133  { \
134  dy_dx[i] = y_sec[i] - x[i] * dy_dv_sec[i] * v / (x[i] * x[i]); \
135  for (unsigned int j = 0; j < _n_secondary_vapors; j++) \
136  { \
137  if (j == i) \
138  continue; \
139  if (_n_secondary_vapors > 1) \
140  imperfectJacobianMessage( \
141  "The mass fraction derivatives in the following " \
142  "function have not been tested for mixtures of 3 or more components:\n\n", \
143  __PRETTY_FUNCTION__); \
144  const Real dxj_dxi = 0; \
145  dy_dx[i] += dxj_dxi * y_sec[j] - x[j] * dy_dv_sec[j] * v / (x[j] * x[j]) * dxj_dxi; \
146  } \
147  const Real dx_primary_dxi = -x_primary / (1. - x[i]); \
148  dy_dx[i] += dx_primary_dxi * y_primary - \
149  x_primary * dy_dv_primary * v / (x_primary * x_primary) * dx_primary_dxi; \
150  } \
151  }
152 
156 #define define_transport_prop_from_T_v(prop) \
157  Real IdealRealGasMixtureFluidProperties::prop##_from_T_v( \
158  Real T, Real v, const std::vector<Real> & x) const \
159  { \
160  const Real x_primary = primaryMassFraction(x); \
161  Real M_primary = _fp_primary->molarMass(); \
162  \
163  Real sum = x_primary / M_primary; \
164  for (unsigned int i = 0; i < _n_secondary_vapors; i++) \
165  sum += x[i] / _fp_secondary[i]->molarMass(); \
166  const Real M_star = 1. / sum; \
167  \
168  const Real vp = v / x_primary; \
169  const Real ep = _fp_primary->e_from_T_v(T, vp); \
170  const Real yp = _fp_primary->prop##_from_v_e(vp, ep); \
171  Real y = x_primary * M_star / M_primary * yp; \
172  \
173  for (unsigned int i = 0; i < _n_secondary_vapors; i++) \
174  { \
175  const Real vi = v / x[i]; \
176  const Real ei = _fp_secondary[i]->e_from_T_v(T, vp); \
177  const Real Mi = _fp_secondary[i]->molarMass(); \
178  const Real yi = _fp_secondary[i]->prop##_from_v_e(vi, ei); \
179  y += x[i] * M_star / Mi * yi; \
180  } \
181  \
182  return y; \
183  } \
184  \
185  void IdealRealGasMixtureFluidProperties::prop##_from_T_v(Real T, \
186  Real v, \
187  const std::vector<Real> & x, \
188  Real & y, \
189  Real & dy_dT, \
190  Real & dy_dv, \
191  std::vector<Real> & dy_dx) const \
192  { \
193  const Real x_primary = primaryMassFraction(x); \
194  Real M_primary = _fp_primary->molarMass(); \
195  \
196  Real sum = x_primary / M_primary; \
197  for (unsigned int i = 0; i < _n_secondary_vapors; i++) \
198  sum += x[i] / _fp_secondary[i]->molarMass(); \
199  const Real M_star = 1. / sum; \
200  \
201  const Real vp = v / x_primary; \
202  const Real ep = _fp_primary->e_from_T_v(T, vp); \
203  const Real yp = _fp_primary->prop##_from_v_e(vp, ep); \
204  y = x_primary * M_star / M_primary * yp; \
205  \
206  imperfectJacobianMessage("The temperature and specific volume derivatives in the following " \
207  "function are currently neglected:\n\n", \
208  __PRETTY_FUNCTION__); \
209  dy_dT = 0; \
210  dy_dv = 0; \
211  \
212  Real sum_yj = 0; \
213  dy_dx.resize(_n_secondary_vapors); \
214  for (unsigned int i = 0; i < _n_secondary_vapors; i++) \
215  { \
216  const Real vi = v / x[i]; \
217  const Real ei = _fp_secondary[i]->e_from_T_v(T, vi); \
218  const Real Mi = _fp_secondary[i]->molarMass(); \
219  const Real yi = _fp_secondary[i]->prop##_from_v_e(vi, ei); \
220  y += x[i] * M_star / Mi * yi; \
221  dy_dx[i] = -M_star / M_primary * yp - \
222  x_primary * (1 / Mi - 1 / M_primary) * M_star * M_star / M_primary * yp + \
223  M_star / Mi * yi; \
224  sum_yj += x[i] * yi / Mi; \
225  } \
226  \
227  for (unsigned int i = 0; i < _n_secondary_vapors; i++) \
228  { \
229  const Real Mi = _fp_secondary[i]->molarMass(); \
230  dy_dx[i] += -M_star * M_star * (1 / Mi - 1 / M_primary) * sum_yj; \
231  } \
232  }
233 
234 // clang-format off
237 
238 define_transport_prop_from_T_v(mu)
239 define_transport_prop_from_T_v(k)
240 
241 define_from_p_T_using_T_v(e)
242 define_from_p_T_using_T_v(s)
243 define_from_p_T_using_T_v(c)
244 define_from_p_T_using_T_v(cp)
245 define_from_p_T_using_T_v(cv)
246 define_from_p_T_using_T_v(mu)
247 define_from_p_T_using_T_v(k)
248 
249 define_from_v_e_using_T_v(p)
250 define_from_v_e_using_T_v(c)
251 
253 // clang-format on
254 {
256  params += NaNInterface::validParams();
257 
258  params.addClassDescription("Class for fluid properties of an arbitrary vapor mixture");
259 
260  params.addRequiredParam<UserObjectName>(
261  "fp_primary", "Name of fluid properties user object for primary vapor component");
262  params.addRequiredParam<std::vector<UserObjectName>>(
263  "fp_secondary", "Name of fluid properties user object(s) for secondary vapor component(s)");
264  params.addParam<Real>("_T_mix_max", 1300., "Maximum temperature of the mixture");
265 
266  return params;
267 }
268 
270  const InputParameters & parameters)
271  : VaporMixtureFluidProperties(parameters),
272  NaNInterface(this),
273  _fp_primary(&getUserObject<SinglePhaseFluidProperties>("fp_primary")),
274  _fp_secondary_names(getParam<std::vector<UserObjectName>>("fp_secondary")),
275  _n_secondary_vapors(_fp_secondary_names.size()),
276  _T_mix_max(getParam<Real>("_T_mix_max"))
277 {
279  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
280  _fp_secondary[i] = &getUserObjectByName<SinglePhaseFluidProperties>(_fp_secondary_names[i]);
281 }
282 
285 {
286  return *_fp_primary;
287 }
288 
291 {
292  mooseAssert(i < getNumberOfSecondaryVapors(), "Requested secondary index too high.");
293  return *_fp_secondary[i];
294 }
295 
296 Real
297 IdealRealGasMixtureFluidProperties::T_from_v_e(Real v, Real e, const std::vector<Real> & x) const
298 {
299  Real v_primary = v / primaryMassFraction(x);
300  static const Real vc = 1. / _fp_primary->criticalDensity();
301  static const Real ec = _fp_primary->criticalInternalEnergy();
302 
303  // Initial estimate of a bracketing interval for the temperature
304  Real lower_temperature, upper_temperature;
305  if (v_primary > vc)
306  {
307  Real e_sat_primary = _fp_primary->e_spndl_from_v(v_primary);
308  lower_temperature = _fp_primary->T_from_v_e(v_primary, e_sat_primary);
309  }
310  else
311  lower_temperature = _fp_primary->T_from_v_e(v_primary, ec);
312 
313  upper_temperature = _T_mix_max;
314 
315  // Use BrentsMethod to find temperature
316  auto energy_diff = [&v, &e, &x, this](Real T) { return this->e_from_T_v(T, v, x) - e; };
317 
318  BrentsMethod::bracket(energy_diff, lower_temperature, upper_temperature);
319  return BrentsMethod::root(energy_diff, lower_temperature, upper_temperature);
320 }
321 
322 void
323 IdealRealGasMixtureFluidProperties::T_from_v_e(Real v,
324  Real e,
325  const std::vector<Real> & x,
326  Real & T,
327  Real & dT_dv,
328  Real & dT_de,
329  std::vector<Real> & dT_dx) const
330 {
331  T = T_from_v_e(v, e, x);
332 
333  Real e_unused, de_dT_v, de_dv_T;
334  std::vector<Real> de_dx_Tv;
335  e_from_T_v(T, v, x, e_unused, de_dT_v, de_dv_T, de_dx_Tv);
336 
337  dT_dv = -de_dv_T / de_dT_v;
338  dT_de = 1.0 / de_dT_v;
339  dT_dx.resize(_n_secondary_vapors);
340  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
341  dT_dx[i] = -de_dx_Tv[i] / de_dT_v;
342 }
343 
344 Real
345 IdealRealGasMixtureFluidProperties::rho_from_p_T(Real p, Real T, const std::vector<Real> & x) const
346 {
347  return 1.0 / v_from_p_T(p, T, x);
348 }
349 
350 void
351 IdealRealGasMixtureFluidProperties::rho_from_p_T(Real p,
352  Real T,
353  const std::vector<Real> & x,
354  Real & rho,
355  Real & drho_dp,
356  Real & drho_dT,
357  std::vector<Real> & drho_dx) const
358 {
359  Real v, dv_dp, dv_dT;
360  std::vector<Real> dv_dx;
361  v_from_p_T(p, T, x, v, dv_dp, dv_dT, dv_dx);
362 
363  rho = 1.0 / v;
364  const Real drho_dv = -1.0 / (v * v);
365  drho_dp = drho_dv * dv_dp;
366  drho_dT = drho_dv * dv_dT;
367  drho_dx.resize(_n_secondary_vapors);
368  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
369  drho_dx[i] = drho_dv * dv_dx[i];
370 }
371 
372 Real
373 IdealRealGasMixtureFluidProperties::v_from_p_T(Real p, Real T, const std::vector<Real> & x) const
374 {
375  const Real x_primary = primaryMassFraction(x);
376  Real M_primary = _fp_primary->molarMass();
377 
378  Real sum = x_primary / M_primary;
379  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
380  sum += x[i] / _fp_secondary[i]->molarMass();
381  Real M_star = 1. / sum;
382  Real v_ideal = R_molar * T / (M_star * p);
383 
384  // check range of validity for primary (condensable) component
385  static const Real Tc = _fp_primary->criticalTemperature();
386  static const Real vc = 1. / _fp_primary->criticalDensity();
387  Real v_spndl, e_spndl;
388  if (T < Tc)
389  _fp_primary->v_e_spndl_from_T(T, v_spndl, e_spndl);
390  else
391  v_spndl = vc;
392 
393  Real lower_spec_volume = v_spndl * x_primary;
394  Real upper_spec_volume = v_ideal; // p*v/(RT) <= 1
395 
396  // Initial estimate of a bracketing interval for the temperature
397  Real p_max = p_from_T_v(T, lower_spec_volume, x);
398  if (p > p_max || upper_spec_volume < lower_spec_volume)
399  return getNaN();
400 
401  // Use BrentsMethod to find temperature
402  auto pressure_diff = [&T, &p, &x, this](Real v) { return this->p_from_T_v(T, v, x) - p; };
403 
404  BrentsMethod::bracket(pressure_diff, lower_spec_volume, upper_spec_volume);
405  return BrentsMethod::root(pressure_diff, lower_spec_volume, upper_spec_volume);
406 }
407 
408 void
409 IdealRealGasMixtureFluidProperties::v_from_p_T(Real p,
410  Real T,
411  const std::vector<Real> & x,
412  Real & v,
413  Real & dv_dp,
414  Real & dv_dT,
415  std::vector<Real> & dv_dx) const
416 {
417  v = v_from_p_T(p, T, x);
418 
419  Real p_unused, dp_dT, dp_dv;
420  std::vector<Real> dp_dx;
421  p_from_T_v(T, v, x, p_unused, dp_dT, dp_dv, dp_dx);
422 
423  dv_dp = 1. / dp_dv;
424  dv_dT = -dp_dT / dp_dv;
425  dv_dx.resize(_n_secondary_vapors);
426  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
427  dv_dx[i] = -dp_dx[i] / dp_dv;
428 }
429 
430 Real
431 IdealRealGasMixtureFluidProperties::e_from_p_rho(Real p,
432  Real rho,
433  const std::vector<Real> & x) const
434 {
435  Real v = 1. / rho;
436  Real T = T_from_p_v(p, v, x);
437 
438  return e_from_T_v(T, v, x);
439 }
440 
441 void
442 IdealRealGasMixtureFluidProperties::e_from_p_rho(Real p,
443  Real rho,
444  const std::vector<Real> & x,
445  Real & e,
446  Real & de_dp,
447  Real & de_drho,
448  std::vector<Real> & de_dx) const
449 {
450  Real v = 1. / rho;
451  Real T = T_from_p_v(p, v, x);
452  Real p_, dp_dT, dp_dv, de_dT, de_dv;
453  std::vector<Real> dp_dx_Tv;
454  std::vector<Real> de_dx_Tv;
455 
456  p_from_T_v(T, v, x, p_, dp_dT, dp_dv, dp_dx_Tv);
457  e_from_T_v(T, v, x, e, de_dT, de_dv, de_dx_Tv);
458 
459  de_dp = de_dT / dp_dT;
460  de_drho = (-v * v) * (de_dv - de_dT * dp_dv / dp_dT);
461 
462  // Derivatives with respect to mass fractions:
463  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
464  {
465  de_dx[i] = de_dx_Tv[i] - de_dT * dp_dx_Tv[i] / dp_dT;
466  }
467 }
468 
469 Real
470 IdealRealGasMixtureFluidProperties::T_from_p_v(Real p, Real v, const std::vector<Real> & x) const
471 {
472  Real v_primary = v / primaryMassFraction(x);
473  static const Real vc = 1. / _fp_primary->criticalDensity();
474  static const Real ec = _fp_primary->criticalInternalEnergy();
475 
476  // Initial estimate of a bracketing interval for the temperature
477  Real lower_temperature, upper_temperature;
478  if (v_primary > vc)
479  {
480  Real e_sat_primary = _fp_primary->e_spndl_from_v(v_primary);
481  lower_temperature = _fp_primary->T_from_v_e(v_primary, e_sat_primary);
482  }
483  else
484  lower_temperature = _fp_primary->T_from_v_e(v_primary, ec);
485 
486  upper_temperature = _T_mix_max;
487 
488  // Use BrentsMethod to find temperature
489  auto pressure_diff = [&p, &v, &x, this](Real T) { return this->p_from_T_v(T, v, x) - p; };
490 
491  BrentsMethod::bracket(pressure_diff, lower_temperature, upper_temperature);
492  return BrentsMethod::root(pressure_diff, lower_temperature, upper_temperature);
493 }
494 
495 void
496 IdealRealGasMixtureFluidProperties::T_from_p_v(Real p,
497  Real v,
498  const std::vector<Real> & x,
499  Real & T,
500  Real & dT_dp,
501  Real & dT_dv,
502  std::vector<Real> & dT_dx) const
503 {
504  T = T_from_p_v(p, v, x);
505 
506  // pressure and derivatives
507  Real p_unused, dp_dT_v, dp_dv_T;
508  std::vector<Real> dp_dx_Tv;
509  p_from_T_v(T, v, x, p_unused, dp_dT_v, dp_dv_T, dp_dx_Tv);
510 
511  // Compute derivatives using the following rules:
512  dT_dp = 1. / dp_dT_v;
513  dT_dv = -dp_dv_T / dp_dT_v;
514 
515  // Derivatives with respect to mass fractions:
516  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
517  dT_dx[i] = -dp_dx_Tv[i] / dp_dT_v;
518 }
519 
520 Real
521 IdealRealGasMixtureFluidProperties::p_from_T_v(Real T, Real v, const std::vector<Real> & x) const
522 {
523  const Real x_primary = primaryMassFraction(x);
524  Real p = _fp_primary->p_from_T_v(T, v / x_primary);
525 
526  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
527  p += _fp_secondary[i]->p_from_T_v(T, v / x[i]);
528 
529  return p;
530 }
531 
532 void
533 IdealRealGasMixtureFluidProperties::p_from_T_v(Real T,
534  Real v,
535  const std::vector<Real> & x,
536  Real & p,
537  Real & dp_dT,
538  Real & dp_dv,
539  std::vector<Real> & dp_dx) const
540 {
541  const Real x_primary = primaryMassFraction(x);
542 
543  Real p_primary, dp_dT_primary, dp_dv_primary;
544  _fp_primary->p_from_T_v(T, v / x_primary, p_primary, dp_dT_primary, dp_dv_primary);
545  p = p_primary;
546  dp_dT = dp_dT_primary;
547  dp_dv = dp_dv_primary / x_primary;
548 
549  // get the partial pressures and their derivatives first
550  Real p_sec, dp_dT_sec, dxj_dxi;
551  std::vector<Real> dp_dv_sec(_n_secondary_vapors);
552  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
553  {
554  _fp_secondary[i]->p_from_T_v(T, v / x[i], p_sec, dp_dT_sec, dp_dv_sec[i]);
555  p += p_sec;
556  dp_dT += dp_dT_sec;
557  dp_dv += dp_dv_sec[i] / x[i];
558  }
559 
560  // get the composition dependent derivatives of the secondary vapors
561  dp_dx.resize(_n_secondary_vapors);
562  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
563  {
564  dp_dx[i] = -dp_dv_sec[i] * v / (x[i] * x[i]);
565  for (unsigned int j = 0; j < _n_secondary_vapors; j++)
566  {
567  if (j == i)
568  continue;
569  // Note: this was previously as follows, but we were unable to understand
570  // why and this is not currently tested (requires 3 or more components in
571  // the mixture):
572  // dxj_dxi = -x[j] / (1. - x[i]);
573  if (_n_secondary_vapors > 1)
575  "The mass fraction derivatives in the following "
576  "function have not been tested for mixtures of 3 or more components:\n\n",
577  __PRETTY_FUNCTION__);
578  dxj_dxi = 0;
579  dp_dx[i] += -dp_dv_sec[j] * v / (x[j] * x[j]) * dxj_dxi;
580  }
581  dp_dx[i] += -dp_dv_primary * v / (x_primary * x_primary) * (-x_primary / (1. - x[i]));
582  }
583 }
584 
585 Real
586 IdealRealGasMixtureFluidProperties::c_from_T_v(Real T, Real v, const std::vector<Real> & x) const
587 {
588  Real p, dp_dT, dp_dv;
589  std::vector<Real> dp_dx;
590  p_from_T_v(T, v, x, p, dp_dT, dp_dv, dp_dx);
591 
592  Real s, ds_dT, ds_dv;
593  std::vector<Real> ds_dx;
594  s_from_T_v(T, v, x, s, ds_dT, ds_dv, ds_dx);
595 
596  const Real dp_dv_s = dp_dv - dp_dT * ds_dv / ds_dT;
597 
598  if (dp_dv_s >= 0)
599  mooseWarning("c_from_T_v(): dp_dv_s = ", dp_dv_s, ". Should be negative.");
600  return v * std::sqrt(-dp_dv_s);
601 }
602 
603 void
604 IdealRealGasMixtureFluidProperties::c_from_T_v(Real T,
605  Real v,
606  const std::vector<Real> & x,
607  Real & c,
608  Real & dc_dT,
609  Real & dc_dv,
610  std::vector<Real> & dc_dx) const
611 {
612  c = c_from_T_v(T, v, x);
613 
614  // For derived properties, we would need higher order derivatives;
615  // therefore, numerical derivatives are used here.
616  const Real dT = 1.e-6;
617  const Real T_perturbed = T + dT;
618  Real c_perturbed = c_from_T_v(T_perturbed, v, x);
619  dc_dT = (c_perturbed - c) / dT;
620 
621  const Real dv = v * 1.e-6;
622  const Real v_perturbed = v + dv;
623  c_perturbed = c_from_T_v(T, v_perturbed, x);
624  dc_dv = (c_perturbed - c) / dv;
625 
626  dc_dx.resize(_n_secondary_vapors);
627  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
628  {
629  std::vector<Real> x_perturbed(x);
630  const Real dx_i = 1e-6;
631  for (unsigned int j = 0; j < _n_secondary_vapors; j++)
632  {
633  if (j != i)
634  x_perturbed[j] =
635  x[j] * (1.0 - (x[i] + dx_i)) / (1.0 - x[i]); // recalculate new mass fractions
636  }
637  x_perturbed[i] += dx_i;
638  c_perturbed = c_from_T_v(T, v, x_perturbed);
639  dc_dx[i] = ((c_perturbed - c) / dx_i);
640  }
641 }
642 
643 Real
644 IdealRealGasMixtureFluidProperties::cp_from_T_v(Real T, Real v, const std::vector<Real> & x) const
645 {
646  const Real x_primary = primaryMassFraction(x);
647 
648  Real p, dp_dT, dp_dv;
649  std::vector<Real> dp_dx;
650  p_from_T_v(T, v, x, p, dp_dT, dp_dv, dp_dx);
651 
652  Real h, dh_dT, dh_dv;
653  _fp_primary->h_from_T_v(T, v / x_primary, h, dh_dT, dh_dv);
654  const Real cp_primary = dh_dT - dh_dv * dp_dT / dp_dv;
655  Real cp = x_primary * cp_primary;
656 
657  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
658  {
659  _fp_secondary[i]->h_from_T_v(T, v / x[i], h, dh_dT, dh_dv);
660  const Real cp_sec = dh_dT - dh_dv * dp_dT / dp_dv;
661  cp += x[i] * cp_sec;
662  }
663 
664  return cp;
665 }
666 
667 void
668 IdealRealGasMixtureFluidProperties::cp_from_T_v(Real T,
669  Real v,
670  const std::vector<Real> & x,
671  Real & cp,
672  Real & dcp_dT,
673  Real & dcp_dv,
674  std::vector<Real> & dcp_dx) const
675 {
676  const Real x_primary = primaryMassFraction(x);
677 
678  Real p, dp_dT, dp_dv;
679  std::vector<Real> dp_dx;
680  p_from_T_v(T, v, x, p, dp_dT, dp_dv, dp_dx);
681 
682  Real h, dh_dT, dh_dv;
683  _fp_primary->h_from_T_v(T, v / x_primary, h, dh_dT, dh_dv);
684  const Real cp_primary = dh_dT - dh_dv * dp_dT / dp_dv;
685  cp = x_primary * cp_primary;
686 
687  // Neglect these for now. These require higher-order derivatives, so a finite
688  // difference should probably be used, as for sound speed.
689  imperfectJacobianMessage("The temperature and specific volume derivatives in the following "
690  "function are currently neglected:\n\n",
691  __PRETTY_FUNCTION__);
692  dcp_dT = 0;
693  dcp_dv = 0;
694 
695  dcp_dx.resize(_n_secondary_vapors);
696  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
697  {
698  _fp_secondary[i]->h_from_T_v(T, v / x[i], h, dh_dT, dh_dv);
699  const Real cp_sec = dh_dT - dh_dv * dp_dT / dp_dv;
700  cp += x[i] * cp_sec;
701  dcp_dx[i] = cp_sec - cp_primary;
702  }
703 }
704 
705 Real
706 IdealRealGasMixtureFluidProperties::cv_from_T_v(Real T, Real v, const std::vector<Real> & x) const
707 {
708  const Real x_primary = primaryMassFraction(x);
709  Real cv = x_primary * _fp_primary->cv_from_T_v(T, v / x_primary);
710 
711  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
712  cv += x[i] * _fp_secondary[i]->cv_from_T_v(T, v / x[i]);
713 
714  return cv;
715 }
716 
717 void
718 IdealRealGasMixtureFluidProperties::cv_from_T_v(Real T,
719  Real v,
720  const std::vector<Real> & x,
721  Real & cv,
722  Real & dcv_dT,
723  Real & dcv_dv,
724  std::vector<Real> & dcv_dx) const
725 {
726  const Real x_primary = primaryMassFraction(x);
727  const Real cv_primary = _fp_primary->cv_from_T_v(T, v / x_primary);
728  cv = x_primary * cv_primary;
729 
730  // Neglect these for now. These require higher-order derivatives, so a finite
731  // difference should probably be used, as for sound speed.
732  imperfectJacobianMessage("The temperature and specific volume derivatives in the following "
733  "function are currently neglected:\n\n",
734  __PRETTY_FUNCTION__);
735  dcv_dT = 0;
736  dcv_dv = 0;
737 
738  dcv_dx.resize(_n_secondary_vapors);
739  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
740  {
741  const Real cv_sec = _fp_secondary[i]->cv_from_T_v(T, v / x[i]);
742  cv += x[i] * cv_sec;
743  dcv_dx[i] = cv_sec - cv_primary;
744  }
745 }
746 
747 Real
749  Real T,
750  const std::vector<Real> & x) const
751 {
753  Real xs;
754  if (T > T_c)
755  // return 1. to indicate that the primary fluid will not condense for
756  // given (p,T)
757  xs = 1.;
758  else
759  {
760  Real pp_sat = _fp_primary->pp_sat_from_p_T(p, T);
761  if (pp_sat < 0.)
762  {
763  // return 1. to indicate that the primary fluid will not condense for
764  // given (p,T)
765  xs = 1.;
766  return xs;
767  }
768  Real v_primary = _fp_primary->v_from_p_T(pp_sat, T);
769  Real pp_sat_secondary = p - pp_sat;
770 
771  Real v_secondary;
772  if (_n_secondary_vapors == 1)
773  v_secondary = _fp_secondary[0]->v_from_p_T(pp_sat_secondary, T);
774  else
775  {
776  std::vector<Real> x_sec(_n_secondary_vapors);
777  const Real x_primary = primaryMassFraction(x);
778  Real sum = 0.;
779  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
780  {
781  x_sec[i] = x[i] / (1. - x_primary);
782  sum += x_sec[i] / _fp_secondary[i]->molarMass();
783  }
784  Real M_star = 1. / sum;
785  v_secondary = R_molar * T / (M_star * pp_sat_secondary);
786  int it = 0;
787  double f = 1., df_dvs, pp_sec, p_sec, dp_dT_sec, dp_dv_sec, dp_dv;
788  double tol_p = 1.e-8;
789  while (std::fabs(f / pp_sat_secondary) > tol_p)
790  {
791  pp_sec = 0.;
792  dp_dv = 0.;
793  for (unsigned int i = 0; i < _n_secondary_vapors; i++)
794  {
795  _fp_secondary[i]->p_from_T_v(T, v_secondary / x_sec[i], p_sec, dp_dT_sec, dp_dv_sec);
796  pp_sec += p_sec;
797  dp_dv += dp_dv_sec / x_sec[i];
798  }
799  f = pp_sec - pp_sat_secondary;
800  df_dvs = dp_dv;
801  v_secondary -= f / df_dvs;
802  if (it++ > 15)
803  return getNaN();
804  }
805  }
806 
807  xs = v_secondary / (v_primary + v_secondary);
808  }
809 
810  return xs;
811 }
virtual const SinglePhaseFluidProperties & getSecondaryFluidProperties(unsigned int i=0) const override
Gets a secondary component single-phase fluid properties.
Base class for fluid properties of vapor mixtures.
static const std::string cv
Definition: NS.h:122
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
virtual Real molarMass() const
Molar mass [kg/mol].
Real xs_prim_from_p_T(Real p, Real T, const std::vector< Real > &x) const
Mass fraction of primary (condensable) component at saturation from pressure and temperature.
Real getNaN() const
Throws an error or returns a NaN with or without a warning, with a default message.
Definition: NaNInterface.h:46
void imperfectJacobianMessage(Args... args) const
const SinglePhaseFluidProperties *const _fp_primary
Primary vapor fluid properties.
void mooseWarning(Args &&... args) const
static InputParameters validParams()
void addRequiredParam(const std::string &name, const std::string &doc_string)
const Real _T_mix_max
maximum temperature of all components
virtual Real criticalInternalEnergy() const
Critical specific internal energy.
virtual const SinglePhaseFluidProperties & getPrimaryFluidProperties() const override
Gets the primary component single-phase fluid properties.
static const std::string cp
Definition: NS.h:121
virtual Real criticalTemperature() const
Critical temperature.
const std::vector< double > x
Real f(Real x)
Test function for Brents method.
static InputParameters validParams()
static const std::string mu
Definition: NS.h:123
Real root(std::function< Real(Real)> const &f, Real x1, Real x2, Real tol=1.0e-12)
Finds the root of a function using Brent&#39;s method.
Definition: BrentsMethod.C:66
Common class for single phase fluid properties.
virtual void v_e_spndl_from_T(Real T, Real &v, Real &e) const
Specific internal energy from temperature and specific volume.
static constexpr const Real R_molar
molar (or universal) gas constant
virtual Real e_spndl_from_v(Real v) const
Specific internal energy from temperature and specific volume.
const std::vector< UserObjectName > _fp_secondary_names
Names of secondary vapor fluid properties.
virtual Real criticalDensity() const
Critical density.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string v
Definition: NS.h:84
Class for fluid properties of an arbitrary vapor mixture.
Real primaryMassFraction(const std::vector< Real > &x) const
Computes the mass fraction of the primary vapor given mass fractions of the secondary vapors...
void addClassDescription(const std::string &doc_string)
const unsigned int _n_secondary_vapors
Number of secondary vapors.
static InputParameters validParams()
Definition: NaNInterface.C:15
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
IdealRealGasMixtureFluidProperties(const InputParameters &parameters)
Interface class for producing errors, warnings, or just quiet NaNs.
Definition: NaNInterface.h:22
std::vector< const SinglePhaseFluidProperties * > _fp_secondary
Secondary vapor fluid properties.
registerMooseObject("FluidPropertiesApp", IdealRealGasMixtureFluidProperties)
static const std::string k
Definition: NS.h:130
void bracket(std::function< Real(Real)> const &f, Real &x1, Real &x2)
Function to bracket a root of a given function.
Definition: BrentsMethod.C:17