Line data Source code
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 :
10 : #include "IdealRealGasMixtureFluidProperties.h"
11 : #include "SinglePhaseFluidProperties.h"
12 : #include "BrentsMethod.h"
13 : #include <numeric>
14 :
15 : registerMooseObject("FluidPropertiesApp", IdealRealGasMixtureFluidProperties);
16 :
17 : /**
18 : * Defines the value and derivative methods from (v,e) for a property y
19 : * using T(v,e) and y(T,v).
20 : */
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 :
52 : /**
53 : * Defines the value and derivative methods from (p,T) for a property y
54 : * using v(p,T) and y(T,v).
55 : */
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 :
87 : /**
88 : * Defines the value and derivative methods from (T,v) for a mass-specific property y
89 : */
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 :
153 : /**
154 : * Defines the value and derivative methods from (T,v) for a transport property y
155 : */
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
235 784 : define_mass_specific_prop_from_T_v(e)
236 391 : define_mass_specific_prop_from_T_v(s)
237 :
238 154 : define_transport_prop_from_T_v(mu)
239 154 : define_transport_prop_from_T_v(k)
240 :
241 28 : define_from_p_T_using_T_v(e)
242 43 : define_from_p_T_using_T_v(s)
243 44 : define_from_p_T_using_T_v(c)
244 44 : define_from_p_T_using_T_v(cp)
245 44 : define_from_p_T_using_T_v(cv)
246 43 : define_from_p_T_using_T_v(mu)
247 43 : define_from_p_T_using_T_v(k)
248 :
249 27 : define_from_v_e_using_T_v(p)
250 11 : define_from_v_e_using_T_v(c)
251 :
252 119 : InputParameters IdealRealGasMixtureFluidProperties::validParams()
253 : // clang-format on
254 : {
255 119 : InputParameters params = VaporMixtureFluidProperties::validParams();
256 119 : params += NaNInterface::validParams();
257 :
258 119 : params.addClassDescription("Class for fluid properties of an arbitrary vapor mixture");
259 :
260 238 : params.addRequiredParam<UserObjectName>(
261 : "fp_primary", "Name of fluid properties user object for primary vapor component");
262 238 : params.addRequiredParam<std::vector<UserObjectName>>(
263 : "fp_secondary", "Name of fluid properties user object(s) for secondary vapor component(s)");
264 238 : params.addParam<Real>("_T_mix_max", 1300., "Maximum temperature of the mixture");
265 :
266 119 : return params;
267 0 : }
268 :
269 67 : IdealRealGasMixtureFluidProperties::IdealRealGasMixtureFluidProperties(
270 67 : const InputParameters & parameters)
271 : : VaporMixtureFluidProperties(parameters),
272 : NaNInterface(this),
273 67 : _fp_primary(&getUserObject<SinglePhaseFluidProperties>("fp_primary")),
274 201 : _fp_secondary_names(getParam<std::vector<UserObjectName>>("fp_secondary")),
275 67 : _n_secondary_vapors(_fp_secondary_names.size()),
276 201 : _T_mix_max(getParam<Real>("_T_mix_max"))
277 : {
278 67 : _fp_secondary.resize(_n_secondary_vapors);
279 156 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
280 89 : _fp_secondary[i] = &getUserObjectByName<SinglePhaseFluidProperties>(_fp_secondary_names[i]);
281 67 : }
282 :
283 : const SinglePhaseFluidProperties &
284 0 : IdealRealGasMixtureFluidProperties::getPrimaryFluidProperties() const
285 : {
286 0 : return *_fp_primary;
287 : }
288 :
289 : const SinglePhaseFluidProperties &
290 0 : IdealRealGasMixtureFluidProperties::getSecondaryFluidProperties(unsigned int i) const
291 : {
292 : mooseAssert(i < getNumberOfSecondaryVapors(), "Requested secondary index too high.");
293 0 : return *_fp_secondary[i];
294 : }
295 :
296 : Real
297 59 : IdealRealGasMixtureFluidProperties::T_from_v_e(Real v, Real e, const std::vector<Real> & x) const
298 : {
299 59 : Real v_primary = v / primaryMassFraction(x);
300 59 : static const Real vc = 1. / _fp_primary->criticalDensity();
301 59 : 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 59 : if (v_primary > vc)
306 : {
307 59 : Real e_sat_primary = _fp_primary->e_spndl_from_v(v_primary);
308 59 : lower_temperature = _fp_primary->T_from_v_e(v_primary, e_sat_primary);
309 : }
310 : else
311 0 : lower_temperature = _fp_primary->T_from_v_e(v_primary, ec);
312 :
313 59 : upper_temperature = _T_mix_max;
314 :
315 : // Use BrentsMethod to find temperature
316 336 : auto energy_diff = [&v, &e, &x, this](Real T) { return this->e_from_T_v(T, v, x) - e; };
317 :
318 59 : BrentsMethod::bracket(energy_diff, lower_temperature, upper_temperature);
319 118 : return BrentsMethod::root(energy_diff, lower_temperature, upper_temperature);
320 : }
321 :
322 : void
323 3 : 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 3 : 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 3 : e_from_T_v(T, v, x, e_unused, de_dT_v, de_dv_T, de_dx_Tv);
336 :
337 3 : dT_dv = -de_dv_T / de_dT_v;
338 3 : dT_de = 1.0 / de_dT_v;
339 3 : dT_dx.resize(_n_secondary_vapors);
340 6 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
341 3 : dT_dx[i] = -de_dx_Tv[i] / de_dT_v;
342 3 : }
343 :
344 : Real
345 27 : IdealRealGasMixtureFluidProperties::rho_from_p_T(Real p, Real T, const std::vector<Real> & x) const
346 : {
347 27 : return 1.0 / v_from_p_T(p, T, x);
348 : }
349 :
350 : void
351 1 : 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 1 : v_from_p_T(p, T, x, v, dv_dp, dv_dT, dv_dx);
362 :
363 1 : rho = 1.0 / v;
364 1 : const Real drho_dv = -1.0 / (v * v);
365 1 : drho_dp = drho_dv * dv_dp;
366 1 : drho_dT = drho_dv * dv_dT;
367 1 : drho_dx.resize(_n_secondary_vapors);
368 2 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
369 1 : drho_dx[i] = drho_dv * dv_dx[i];
370 1 : }
371 :
372 : Real
373 329 : IdealRealGasMixtureFluidProperties::v_from_p_T(Real p, Real T, const std::vector<Real> & x) const
374 : {
375 329 : const Real x_primary = primaryMassFraction(x);
376 329 : Real M_primary = _fp_primary->molarMass();
377 :
378 329 : Real sum = x_primary / M_primary;
379 676 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
380 347 : sum += x[i] / _fp_secondary[i]->molarMass();
381 329 : Real M_star = 1. / sum;
382 329 : Real v_ideal = R_molar * T / (M_star * p);
383 :
384 : // check range of validity for primary (condensable) component
385 329 : static const Real Tc = _fp_primary->criticalTemperature();
386 329 : static const Real vc = 1. / _fp_primary->criticalDensity();
387 : Real v_spndl, e_spndl;
388 329 : if (T < Tc)
389 87 : _fp_primary->v_e_spndl_from_T(T, v_spndl, e_spndl);
390 : else
391 242 : v_spndl = vc;
392 :
393 329 : Real lower_spec_volume = v_spndl * x_primary;
394 329 : Real upper_spec_volume = v_ideal; // p*v/(RT) <= 1
395 :
396 : // Initial estimate of a bracketing interval for the temperature
397 329 : Real p_max = p_from_T_v(T, lower_spec_volume, x);
398 329 : if (p > p_max || upper_spec_volume < lower_spec_volume)
399 0 : return getNaN();
400 :
401 : // Use BrentsMethod to find temperature
402 4479 : auto pressure_diff = [&T, &p, &x, this](Real v) { return this->p_from_T_v(T, v, x) - p; };
403 :
404 329 : BrentsMethod::bracket(pressure_diff, lower_spec_volume, upper_spec_volume);
405 658 : return BrentsMethod::root(pressure_diff, lower_spec_volume, upper_spec_volume);
406 : }
407 :
408 : void
409 9 : 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 9 : v = v_from_p_T(p, T, x);
418 :
419 : Real p_unused, dp_dT, dp_dv;
420 : std::vector<Real> dp_dx;
421 9 : p_from_T_v(T, v, x, p_unused, dp_dT, dp_dv, dp_dx);
422 :
423 9 : dv_dp = 1. / dp_dv;
424 9 : dv_dT = -dp_dT / dp_dv;
425 9 : dv_dx.resize(_n_secondary_vapors);
426 18 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
427 9 : dv_dx[i] = -dp_dx[i] / dp_dv;
428 9 : }
429 :
430 : Real
431 8 : IdealRealGasMixtureFluidProperties::e_from_p_rho(Real p,
432 : Real rho,
433 : const std::vector<Real> & x) const
434 : {
435 8 : Real v = 1. / rho;
436 8 : Real T = T_from_p_v(p, v, x);
437 :
438 8 : return e_from_T_v(T, v, x);
439 : }
440 :
441 : void
442 1 : 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 1 : Real v = 1. / rho;
451 1 : 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 1 : p_from_T_v(T, v, x, p_, dp_dT, dp_dv, dp_dx_Tv);
457 1 : e_from_T_v(T, v, x, e, de_dT, de_dv, de_dx_Tv);
458 :
459 1 : de_dp = de_dT / dp_dT;
460 1 : de_drho = (-v * v) * (de_dv - de_dT * dp_dv / dp_dT);
461 :
462 : // Derivatives with respect to mass fractions:
463 2 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
464 : {
465 1 : de_dx[i] = de_dx_Tv[i] - de_dT * dp_dx_Tv[i] / dp_dT;
466 : }
467 1 : }
468 :
469 : Real
470 18 : IdealRealGasMixtureFluidProperties::T_from_p_v(Real p, Real v, const std::vector<Real> & x) const
471 : {
472 18 : Real v_primary = v / primaryMassFraction(x);
473 18 : static const Real vc = 1. / _fp_primary->criticalDensity();
474 18 : 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 18 : if (v_primary > vc)
479 : {
480 18 : Real e_sat_primary = _fp_primary->e_spndl_from_v(v_primary);
481 18 : lower_temperature = _fp_primary->T_from_v_e(v_primary, e_sat_primary);
482 : }
483 : else
484 0 : lower_temperature = _fp_primary->T_from_v_e(v_primary, ec);
485 :
486 18 : upper_temperature = _T_mix_max;
487 :
488 : // Use BrentsMethod to find temperature
489 99 : auto pressure_diff = [&p, &v, &x, this](Real T) { return this->p_from_T_v(T, v, x) - p; };
490 :
491 18 : BrentsMethod::bracket(pressure_diff, lower_temperature, upper_temperature);
492 36 : return BrentsMethod::root(pressure_diff, lower_temperature, upper_temperature);
493 : }
494 :
495 : void
496 1 : 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 1 : 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 1 : 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 1 : dT_dp = 1. / dp_dT_v;
513 1 : dT_dv = -dp_dv_T / dp_dT_v;
514 :
515 : // Derivatives with respect to mass fractions:
516 2 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
517 1 : dT_dx[i] = -dp_dx_Tv[i] / dp_dT_v;
518 1 : }
519 :
520 : Real
521 4939 : IdealRealGasMixtureFluidProperties::p_from_T_v(Real T, Real v, const std::vector<Real> & x) const
522 : {
523 4939 : const Real x_primary = primaryMassFraction(x);
524 4939 : Real p = _fp_primary->p_from_T_v(T, v / x_primary);
525 :
526 10058 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
527 5119 : p += _fp_secondary[i]->p_from_T_v(T, v / x[i]);
528 :
529 4939 : return p;
530 : }
531 :
532 : void
533 121 : 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 121 : const Real x_primary = primaryMassFraction(x);
542 :
543 : Real p_primary, dp_dT_primary, dp_dv_primary;
544 121 : _fp_primary->p_from_T_v(T, v / x_primary, p_primary, dp_dT_primary, dp_dv_primary);
545 121 : p = p_primary;
546 121 : dp_dT = dp_dT_primary;
547 121 : 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 121 : std::vector<Real> dp_dv_sec(_n_secondary_vapors);
552 242 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
553 : {
554 121 : _fp_secondary[i]->p_from_T_v(T, v / x[i], p_sec, dp_dT_sec, dp_dv_sec[i]);
555 121 : p += p_sec;
556 121 : dp_dT += dp_dT_sec;
557 121 : dp_dv += dp_dv_sec[i] / x[i];
558 : }
559 :
560 : // get the composition dependent derivatives of the secondary vapors
561 121 : dp_dx.resize(_n_secondary_vapors);
562 242 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
563 : {
564 121 : dp_dx[i] = -dp_dv_sec[i] * v / (x[i] * x[i]);
565 242 : for (unsigned int j = 0; j < _n_secondary_vapors; j++)
566 : {
567 121 : if (j == i)
568 121 : 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 0 : if (_n_secondary_vapors > 1)
574 0 : imperfectJacobianMessage(
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 0 : dp_dx[i] += -dp_dv_sec[j] * v / (x[j] * x[j]) * dxj_dxi;
580 : }
581 121 : dp_dx[i] += -dp_dv_primary * v / (x_primary * x_primary) * (-x_primary / (1. - x[i]));
582 : }
583 121 : }
584 :
585 : Real
586 57 : 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 57 : 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 57 : s_from_T_v(T, v, x, s, ds_dT, ds_dv, ds_dx);
595 :
596 57 : const Real dp_dv_s = dp_dv - dp_dT * ds_dv / ds_dT;
597 :
598 57 : if (dp_dv_s >= 0)
599 0 : mooseWarning("c_from_T_v(): dp_dv_s = ", dp_dv_s, ". Should be negative.");
600 57 : return v * std::sqrt(-dp_dv_s);
601 57 : }
602 :
603 : void
604 2 : 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 2 : 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 2 : const Real T_perturbed = T + dT;
618 2 : Real c_perturbed = c_from_T_v(T_perturbed, v, x);
619 2 : dc_dT = (c_perturbed - c) / dT;
620 :
621 2 : const Real dv = v * 1.e-6;
622 2 : const Real v_perturbed = v + dv;
623 2 : c_perturbed = c_from_T_v(T, v_perturbed, x);
624 2 : dc_dv = (c_perturbed - c) / dv;
625 :
626 2 : dc_dx.resize(_n_secondary_vapors);
627 4 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
628 : {
629 2 : std::vector<Real> x_perturbed(x);
630 : const Real dx_i = 1e-6;
631 4 : for (unsigned int j = 0; j < _n_secondary_vapors; j++)
632 : {
633 2 : if (j != i)
634 0 : x_perturbed[j] =
635 0 : x[j] * (1.0 - (x[i] + dx_i)) / (1.0 - x[i]); // recalculate new mass fractions
636 : }
637 2 : x_perturbed[i] += dx_i;
638 2 : c_perturbed = c_from_T_v(T, v, x_perturbed);
639 2 : dc_dx[i] = ((c_perturbed - c) / dx_i);
640 2 : }
641 2 : }
642 :
643 : Real
644 49 : IdealRealGasMixtureFluidProperties::cp_from_T_v(Real T, Real v, const std::vector<Real> & x) const
645 : {
646 49 : const Real x_primary = primaryMassFraction(x);
647 :
648 : Real p, dp_dT, dp_dv;
649 : std::vector<Real> dp_dx;
650 49 : p_from_T_v(T, v, x, p, dp_dT, dp_dv, dp_dx);
651 :
652 : Real h, dh_dT, dh_dv;
653 49 : _fp_primary->h_from_T_v(T, v / x_primary, h, dh_dT, dh_dv);
654 49 : const Real cp_primary = dh_dT - dh_dv * dp_dT / dp_dv;
655 49 : Real cp = x_primary * cp_primary;
656 :
657 98 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
658 : {
659 49 : _fp_secondary[i]->h_from_T_v(T, v / x[i], h, dh_dT, dh_dv);
660 49 : const Real cp_sec = dh_dT - dh_dv * dp_dT / dp_dv;
661 49 : cp += x[i] * cp_sec;
662 : }
663 :
664 49 : return cp;
665 49 : }
666 :
667 : void
668 2 : 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 2 : const Real x_primary = primaryMassFraction(x);
677 :
678 : Real p, dp_dT, dp_dv;
679 : std::vector<Real> dp_dx;
680 2 : p_from_T_v(T, v, x, p, dp_dT, dp_dv, dp_dx);
681 :
682 : Real h, dh_dT, dh_dv;
683 2 : _fp_primary->h_from_T_v(T, v / x_primary, h, dh_dT, dh_dv);
684 2 : const Real cp_primary = dh_dT - dh_dv * dp_dT / dp_dv;
685 2 : 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 2 : imperfectJacobianMessage("The temperature and specific volume derivatives in the following "
690 : "function are currently neglected:\n\n",
691 : __PRETTY_FUNCTION__);
692 2 : dcp_dT = 0;
693 2 : dcp_dv = 0;
694 :
695 2 : dcp_dx.resize(_n_secondary_vapors);
696 4 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
697 : {
698 2 : _fp_secondary[i]->h_from_T_v(T, v / x[i], h, dh_dT, dh_dv);
699 2 : const Real cp_sec = dh_dT - dh_dv * dp_dT / dp_dv;
700 2 : cp += x[i] * cp_sec;
701 2 : dcp_dx[i] = cp_sec - cp_primary;
702 : }
703 2 : }
704 :
705 : Real
706 49 : IdealRealGasMixtureFluidProperties::cv_from_T_v(Real T, Real v, const std::vector<Real> & x) const
707 : {
708 49 : const Real x_primary = primaryMassFraction(x);
709 49 : Real cv = x_primary * _fp_primary->cv_from_T_v(T, v / x_primary);
710 :
711 98 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
712 49 : cv += x[i] * _fp_secondary[i]->cv_from_T_v(T, v / x[i]);
713 :
714 49 : return cv;
715 : }
716 :
717 : void
718 2 : 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 2 : const Real x_primary = primaryMassFraction(x);
727 2 : const Real cv_primary = _fp_primary->cv_from_T_v(T, v / x_primary);
728 2 : 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 2 : imperfectJacobianMessage("The temperature and specific volume derivatives in the following "
733 : "function are currently neglected:\n\n",
734 : __PRETTY_FUNCTION__);
735 2 : dcv_dT = 0;
736 2 : dcv_dv = 0;
737 :
738 2 : dcv_dx.resize(_n_secondary_vapors);
739 4 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
740 : {
741 2 : const Real cv_sec = _fp_secondary[i]->cv_from_T_v(T, v / x[i]);
742 2 : cv += x[i] * cv_sec;
743 2 : dcv_dx[i] = cv_sec - cv_primary;
744 : }
745 2 : }
746 :
747 : Real
748 0 : IdealRealGasMixtureFluidProperties::xs_prim_from_p_T(Real p,
749 : Real T,
750 : const std::vector<Real> & x) const
751 : {
752 0 : Real T_c = _fp_primary->criticalTemperature();
753 : Real xs;
754 0 : 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 0 : Real pp_sat = _fp_primary->pp_sat_from_p_T(p, T);
761 0 : 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 0 : Real v_primary = _fp_primary->v_from_p_T(pp_sat, T);
769 0 : Real pp_sat_secondary = p - pp_sat;
770 :
771 : Real v_secondary;
772 0 : if (_n_secondary_vapors == 1)
773 0 : v_secondary = _fp_secondary[0]->v_from_p_T(pp_sat_secondary, T);
774 : else
775 : {
776 0 : std::vector<Real> x_sec(_n_secondary_vapors);
777 0 : const Real x_primary = primaryMassFraction(x);
778 : Real sum = 0.;
779 0 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
780 : {
781 0 : x_sec[i] = x[i] / (1. - x_primary);
782 0 : sum += x_sec[i] / _fp_secondary[i]->molarMass();
783 : }
784 0 : Real M_star = 1. / sum;
785 0 : 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 0 : while (std::fabs(f / pp_sat_secondary) > tol_p)
790 : {
791 : pp_sec = 0.;
792 : dp_dv = 0.;
793 0 : for (unsigned int i = 0; i < _n_secondary_vapors; i++)
794 : {
795 0 : _fp_secondary[i]->p_from_T_v(T, v_secondary / x_sec[i], p_sec, dp_dT_sec, dp_dv_sec);
796 0 : pp_sec += p_sec;
797 0 : dp_dv += dp_dv_sec / x_sec[i];
798 : }
799 0 : f = pp_sec - pp_sat_secondary;
800 : df_dvs = dp_dv;
801 0 : v_secondary -= f / df_dvs;
802 0 : if (it++ > 15)
803 0 : return getNaN();
804 : }
805 0 : }
806 :
807 0 : xs = v_secondary / (v_primary + v_secondary);
808 : }
809 :
810 : return xs;
811 : }
|
