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 "StiffenedGasTwoPhaseFluidProperties.h"
11 : #include "StiffenedGasFluidProperties.h"
12 :
13 : registerMooseObject("FluidPropertiesApp", StiffenedGasTwoPhaseFluidProperties);
14 :
15 : InputParameters
16 88 : StiffenedGasTwoPhaseFluidProperties::validParams()
17 : {
18 88 : InputParameters params = TwoPhaseFluidProperties::validParams();
19 88 : params += NaNInterface::validParams();
20 :
21 176 : params.addParam<Real>("gamma_liquid", 2.35, "Liquid heat capacity ratio");
22 176 : params.addParam<Real>("cv_liquid", 1816.0, "Liquid isochoric specific heat capacity");
23 176 : params.addParam<Real>("q_liquid", -1.167e6, "Liquid reference specific internal energy");
24 176 : params.addParam<Real>("p_inf_liquid", 1.0e9, "Liquid stiffness pressure");
25 176 : params.addParam<Real>("q_prime_liquid", 0, "Liquid reference specific entropy");
26 176 : params.addParam<Real>("k_liquid", 0.5, "Liquid thermal conductivity");
27 176 : params.addParam<Real>("mu_liquid", 281.8e-6, "Liquid dynamic viscosity");
28 176 : params.addParam<Real>("M_liquid", 0.01801488, "Liquid molar mass");
29 :
30 176 : params.addParam<Real>("gamma_vapor", 1.43, "Vapor heat capacity ratio");
31 176 : params.addParam<Real>("cv_vapor", 1040.0, "Vapor isochoric specific heat capacity");
32 176 : params.addParam<Real>("q_vapor", 2.03e6, "Vapor reference specific internal energy");
33 176 : params.addParam<Real>("p_inf_vapor", 0.0, "Vapor stiffness pressure");
34 176 : params.addParam<Real>("q_prime_vapor", -2.3e4, "Vapor reference specific entropy");
35 176 : params.addParam<Real>("k_vapor", 0.026, "Vapor thermal conductivity");
36 176 : params.addParam<Real>("mu_vapor", 134.4e-7, "Vapor dynamic viscosity");
37 176 : params.addParam<Real>("M_vapor", 0.01801488, "Vapor molar mass");
38 :
39 176 : params.addParam<Real>("T_c", 647.096, "Critical temperature [K]");
40 176 : params.addParam<Real>("p_c", 22.09e6, "Critical pressure [Pa]");
41 176 : params.addParam<Real>("rho_c", 322.0, "Critical density [kg/m^3]");
42 176 : params.addParam<Real>("e_c", 2702979.84310559, "Critical specific internal energy [J/kg]");
43 176 : params.addParam<Real>("T_triple", 273.16, "Triple-point temperature [K]");
44 176 : params.addParam<Real>("L_fusion", 0.334, "Latent heat of fusion [J/kg]");
45 :
46 176 : params.addParam<Real>(
47 176 : "sigma_A", 0.2358, "'A' constant used in surface tension correlation [N/m]");
48 176 : params.addParam<Real>("sigma_B", 1.256, "'B' constant used in surface tension correlation");
49 176 : params.addParam<Real>("sigma_C", 0.625, "'C' constant used in surface tension correlation");
50 :
51 : // Default values correspond to water, from the freezing point to critical point, with 1 K
52 : // increments
53 176 : params.addParam<Real>("T_sat_min", 274.0, "Minimum temperature value in saturation curve [K]");
54 176 : params.addParam<Real>("T_sat_max", 647.0, "Maximum temperature value in saturation curve [K]");
55 176 : params.addParam<Real>(
56 176 : "p_sat_guess", 611.0, "Initial guess for saturation pressure Newton solve [Pa]");
57 176 : params.addParam<unsigned int>(
58 176 : "n_sat_samples", 374, "Number of samples to take in saturation curve");
59 :
60 88 : params.addClassDescription("Two-phase stiffened gas fluid properties");
61 :
62 88 : return params;
63 0 : }
64 :
65 50 : StiffenedGasTwoPhaseFluidProperties::StiffenedGasTwoPhaseFluidProperties(
66 50 : const InputParameters & parameters)
67 : : TwoPhaseFluidProperties(parameters),
68 : NaNInterface(this),
69 :
70 50 : _gamma_liquid(getParam<Real>("gamma_liquid")),
71 100 : _cv_liquid(getParam<Real>("cv_liquid")),
72 50 : _cp_liquid(_gamma_liquid * _cv_liquid),
73 100 : _q_liquid(getParam<Real>("q_liquid")),
74 100 : _p_inf_liquid(getParam<Real>("p_inf_liquid")),
75 100 : _q_prime_liquid(getParam<Real>("q_prime_liquid")),
76 :
77 100 : _gamma_vapor(getParam<Real>("gamma_vapor")),
78 100 : _cv_vapor(getParam<Real>("cv_vapor")),
79 50 : _cp_vapor(_gamma_vapor * _cv_vapor),
80 100 : _q_vapor(getParam<Real>("q_vapor")),
81 100 : _p_inf_vapor(getParam<Real>("p_inf_vapor")),
82 100 : _q_prime_vapor(getParam<Real>("q_prime_vapor")),
83 :
84 100 : _T_c(getParam<Real>("T_c")),
85 100 : _p_c(getParam<Real>("p_c")),
86 100 : _T_triple(getParam<Real>("T_triple")),
87 100 : _L_fusion(getParam<Real>("L_fusion")),
88 :
89 100 : _sigma_A(getParam<Real>("sigma_A")),
90 100 : _sigma_B(getParam<Real>("sigma_B")),
91 100 : _sigma_C(getParam<Real>("sigma_C")),
92 :
93 100 : _T_sat_min(getParam<Real>("T_sat_min")),
94 100 : _T_sat_max(getParam<Real>("T_sat_max")),
95 100 : _p_sat_guess(getParam<Real>("p_sat_guess")),
96 100 : _n_sat_samples(getParam<unsigned int>("n_sat_samples")),
97 50 : _dT_sat((_T_sat_max - _T_sat_min) / (_n_sat_samples - 1)),
98 :
99 50 : _A((_cp_liquid - _cp_vapor + _q_prime_vapor - _q_prime_liquid) / (_cp_vapor - _cv_vapor)),
100 50 : _B((_q_liquid - _q_vapor) / (_cp_vapor - _cv_vapor)),
101 50 : _C((_cp_vapor - _cp_liquid) / (_cp_vapor - _cv_vapor)),
102 50 : _D((_cp_liquid - _cv_liquid) / (_cp_vapor - _cv_vapor)),
103 :
104 50 : _newton_tol(1.0e-8),
105 100 : _newton_max_iter(200)
106 : {
107 50 : if (_tid == 0)
108 : {
109 38 : std::string class_name = "StiffenedGasFluidProperties";
110 38 : InputParameters params = _app.getFactory().getValidParams(class_name);
111 114 : params.set<MooseEnum>("emit_on_nan") = getParam<MooseEnum>("emit_on_nan");
112 38 : params.set<bool>("allow_nonphysical_states") = false;
113 38 : params.set<Real>("gamma") = _gamma_liquid;
114 38 : params.set<Real>("cv") = _cv_liquid;
115 38 : params.set<Real>("q") = _q_liquid;
116 38 : params.set<Real>("p_inf") = _p_inf_liquid;
117 38 : params.set<Real>("q_prime") = _q_prime_liquid;
118 76 : params.set<Real>("k") = getParam<Real>("k_liquid");
119 76 : params.set<Real>("mu") = getParam<Real>("mu_liquid");
120 76 : params.set<Real>("M") = getParam<Real>("M_liquid");
121 38 : _fe_problem.addUserObject(class_name, _liquid_name, params);
122 38 : }
123 50 : _fp_liquid = &_fe_problem.getUserObject<SinglePhaseFluidProperties>(_liquid_name, _tid);
124 :
125 50 : if (_tid == 0)
126 : {
127 38 : std::string class_name = "StiffenedGasFluidProperties";
128 38 : InputParameters params = _app.getFactory().getValidParams(class_name);
129 114 : params.set<MooseEnum>("emit_on_nan") = getParam<MooseEnum>("emit_on_nan");
130 38 : params.set<bool>("allow_nonphysical_states") = false;
131 38 : params.set<Real>("gamma") = _gamma_vapor;
132 38 : params.set<Real>("cv") = _cv_vapor;
133 38 : params.set<Real>("q") = _q_vapor;
134 38 : params.set<Real>("p_inf") = _p_inf_vapor;
135 38 : params.set<Real>("q_prime") = _q_prime_vapor;
136 76 : params.set<Real>("k") = getParam<Real>("k_vapor");
137 76 : params.set<Real>("mu") = getParam<Real>("mu_vapor");
138 76 : params.set<Real>("M") = getParam<Real>("M_vapor");
139 76 : params.set<Real>("T_c") = getParam<Real>("T_c");
140 76 : params.set<Real>("rho_c") = getParam<Real>("rho_c");
141 76 : params.set<Real>("e_c") = getParam<Real>("e_c");
142 38 : _fe_problem.addUserObject(class_name, _vapor_name, params);
143 38 : }
144 50 : _fp_vapor = &_fe_problem.getUserObject<SinglePhaseFluidProperties>(_vapor_name, _tid);
145 :
146 : // calculate the saturation curve p(T) and store the data in two vectors for re-use
147 : {
148 50 : _T_vec.resize(_n_sat_samples);
149 50 : _p_sat_vec.resize(_n_sat_samples);
150 :
151 : // loop over sample temperatures, starting with the minimum
152 50 : Real T = _T_sat_min;
153 18750 : for (unsigned int i = 0; i < _n_sat_samples; i++)
154 : {
155 18700 : _T_vec[i] = T;
156 18700 : _p_sat_vec[i] = compute_p_sat(T);
157 :
158 : // increment sample temperature
159 18700 : T += _dT_sat;
160 : }
161 : }
162 :
163 50 : _ipol_pressure.setData(_T_vec, _p_sat_vec);
164 50 : _ipol_temp.setData(_p_sat_vec, _T_vec);
165 50 : }
166 :
167 : Real
168 1 : StiffenedGasTwoPhaseFluidProperties::p_critical() const
169 : {
170 1 : return _p_c;
171 : }
172 :
173 : Real
174 1 : StiffenedGasTwoPhaseFluidProperties::T_triple() const
175 : {
176 1 : return _T_triple;
177 : }
178 :
179 : Real
180 1 : StiffenedGasTwoPhaseFluidProperties::L_fusion() const
181 : {
182 1 : return _L_fusion;
183 : }
184 :
185 : Real
186 11 : StiffenedGasTwoPhaseFluidProperties::T_sat(Real pressure) const
187 : {
188 11 : return _ipol_temp.sample(pressure);
189 : }
190 :
191 : Real
192 15 : StiffenedGasTwoPhaseFluidProperties::p_sat(Real temperature) const
193 : {
194 15 : return _ipol_pressure.sample(temperature);
195 : }
196 :
197 : Real
198 18700 : StiffenedGasTwoPhaseFluidProperties::compute_p_sat(const Real & T) const
199 : {
200 18700 : Real p_sat = _p_sat_guess;
201 : Real f_norm = 1.0e5;
202 : unsigned int iter = 1;
203 199400 : while (std::fabs(f_norm / p_sat) > _newton_tol)
204 : {
205 : // check for maximum iteration
206 180700 : if (iter > _newton_max_iter)
207 0 : mooseError("StiffenedGasTwoPhaseFluidProperties::compute_p_sat: ",
208 : "Newton solve did not converge after ",
209 0 : _newton_max_iter,
210 : " iterations.");
211 :
212 : // evaluate nonlinear residual and its derivative w.r.t. p_sat
213 180700 : Real f = std::log(p_sat + _p_inf_vapor) - _A - _B / T - _C * std::log(T) -
214 180700 : _D * std::log(p_sat + _p_inf_liquid);
215 180700 : Real f_prime = 1.0 / (p_sat + _p_inf_vapor) - _D / (p_sat + _p_inf_liquid);
216 :
217 : // take Newton step
218 180700 : f_norm = f / f_prime;
219 180700 : p_sat -= f_norm;
220 :
221 180700 : iter++;
222 : }
223 18700 : return p_sat;
224 : }
225 :
226 : Real
227 17 : StiffenedGasTwoPhaseFluidProperties::dT_sat_dp(Real pressure) const
228 : {
229 17 : return _ipol_temp.sampleDerivative(pressure);
230 : }
231 :
232 : Real
233 11 : StiffenedGasTwoPhaseFluidProperties::sigma_from_T(Real T) const
234 : {
235 11 : const Real aux = 1.0 - T / _T_c;
236 11 : return _sigma_A * std::pow(aux, _sigma_B) * (1.0 - _sigma_C * aux);
237 : }
238 :
239 : Real
240 9 : StiffenedGasTwoPhaseFluidProperties::dsigma_dT_from_T(Real T) const
241 : {
242 9 : const Real aux = 1.0 - T / _T_c;
243 9 : const Real daux_dT = -1.0 / _T_c;
244 : const Real dsigma_daux =
245 9 : _sigma_A * (_sigma_B * std::pow(aux, _sigma_B - 1.0) * (1.0 - _sigma_C * aux) -
246 9 : _sigma_C * std::pow(aux, _sigma_B));
247 9 : return dsigma_daux * daux_dT;
248 : }
|
