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