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 "ADShaftConnectedCompressor1PhaseUserObject.h"
11 : #include "SinglePhaseFluidProperties.h"
12 : #include "VolumeJunction1Phase.h"
13 : #include "MooseVariableScalar.h"
14 : #include "THMIndicesVACE.h"
15 : #include "Function.h"
16 : #include "Numerics.h"
17 : #include "metaphysicl/parallel_numberarray.h"
18 : #include "metaphysicl/parallel_dualnumber.h"
19 : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
20 : #include "libmesh/parallel_algebra.h"
21 :
22 : registerMooseObject("ThermalHydraulicsApp", ADShaftConnectedCompressor1PhaseUserObject);
23 :
24 : InputParameters
25 117 : ADShaftConnectedCompressor1PhaseUserObject::validParams()
26 : {
27 117 : InputParameters params = ADVolumeJunction1PhaseUserObject::validParams();
28 117 : params += ADShaftConnectableUserObjectInterface::validParams();
29 :
30 234 : params.addParam<BoundaryName>("inlet", "Compressor inlet");
31 234 : params.addParam<BoundaryName>("outlet", "Compressor outlet");
32 234 : params.addRequiredParam<Point>("di_out", "Direction of connected outlet");
33 234 : params.addRequiredParam<bool>("treat_as_turbine", "Treat the compressor as a turbine?");
34 234 : params.addRequiredParam<Real>("omega_rated", "Rated compressor speed [rad/s]");
35 234 : params.addRequiredParam<Real>("mdot_rated", "Rated compressor mass flow rate [kg/s]");
36 234 : params.addRequiredParam<Real>("rho0_rated",
37 : "Rated compressor inlet stagnation fluid density [kg/m^3]");
38 234 : params.addRequiredParam<Real>("c0_rated", "Rated compressor inlet stagnation sound speed [m/s]");
39 234 : params.addRequiredParam<Real>("speed_cr_fr", "Compressor speed threshold for friction [-]");
40 234 : params.addRequiredParam<Real>("tau_fr_const", "Compressor friction constant [N-m]");
41 234 : params.addRequiredParam<std::vector<Real>>("tau_fr_coeff", "Friction coefficients [N-m]");
42 234 : params.addRequiredParam<Real>("speed_cr_I", "Compressor speed threshold for inertia [-]");
43 234 : params.addRequiredParam<Real>("inertia_const", "Compressor inertia constant [kg-m^2]");
44 234 : params.addRequiredParam<std::vector<Real>>("inertia_coeff",
45 : "Compressor inertia coefficients [kg-m^2]");
46 234 : params.addRequiredParam<std::vector<Real>>(
47 : "speeds",
48 : "Relative corrected speeds. Order of speeds needs correspond to the "
49 : "orders of `Rp_functions` and `eff_functions` [-]");
50 234 : params.addRequiredParam<std::vector<FunctionName>>(
51 : "Rp_functions",
52 : "Functions of pressure ratio versus relative corrected flow. Each function is for a "
53 : "different, constant relative corrected speed. The order of function names should correspond "
54 : "to the order of speeds in the `speeds` parameter [-]");
55 234 : params.addRequiredParam<std::vector<FunctionName>>(
56 : "eff_functions",
57 : "Functions of adiabatic efficiency versus relative corrected flow. Each function is for a "
58 : "different, constant relative corrected speed. The order of function names should correspond "
59 : "to the order of speeds in the `speeds` parameter [-]");
60 234 : params.addRequiredParam<Real>("min_pressure_ratio", "Minimum pressure ratio");
61 234 : params.addRequiredParam<Real>("max_pressure_ratio", "Maximum pressure ratio");
62 234 : params.addRequiredParam<std::string>("compressor_name",
63 : "Name of the instance of this compressor component");
64 234 : params.addRequiredCoupledVar("omega", "Shaft rotational speed [rad/s]");
65 :
66 117 : params.addClassDescription("Computes and caches flux and residual vectors for a 1-phase "
67 : "compressor. Also computes compressor torque "
68 : "and delta_p which is passed to the connected shaft");
69 :
70 117 : return params;
71 0 : }
72 :
73 62 : ADShaftConnectedCompressor1PhaseUserObject::ADShaftConnectedCompressor1PhaseUserObject(
74 62 : const InputParameters & params)
75 : : ADVolumeJunction1PhaseUserObject(params),
76 : ADShaftConnectableUserObjectInterface(this),
77 :
78 62 : _di_out(getParam<Point>("di_out")),
79 124 : _treat_as_turbine(getParam<bool>("treat_as_turbine")),
80 124 : _omega_rated(getParam<Real>("omega_rated")),
81 124 : _mdot_rated(getParam<Real>("mdot_rated")),
82 124 : _rho0_rated(getParam<Real>("rho0_rated")),
83 124 : _c0_rated(getParam<Real>("c0_rated")),
84 124 : _speed_cr_fr(getParam<Real>("speed_cr_fr")),
85 124 : _tau_fr_const(getParam<Real>("tau_fr_const")),
86 124 : _tau_fr_coeff(getParam<std::vector<Real>>("tau_fr_coeff")),
87 124 : _speed_cr_I(getParam<Real>("speed_cr_I")),
88 124 : _inertia_const(getParam<Real>("inertia_const")),
89 124 : _inertia_coeff(getParam<std::vector<Real>>("inertia_coeff")),
90 124 : _speeds(getParam<std::vector<Real>>("speeds")),
91 124 : _Rp_function_names(getParam<std::vector<FunctionName>>("Rp_functions")),
92 124 : _eff_function_names(getParam<std::vector<FunctionName>>("eff_functions")),
93 62 : _n_speeds(_speeds.size()),
94 62 : _Rp_functions(_n_speeds),
95 62 : _eff_functions(_n_speeds),
96 124 : _Rp_min(getParam<Real>("min_pressure_ratio")),
97 124 : _Rp_max(getParam<Real>("max_pressure_ratio")),
98 124 : _compressor_name(getParam<std::string>("compressor_name")),
99 124 : _omega(adCoupledScalarValue("omega"))
100 : {
101 62 : if (_n_speeds != _Rp_function_names.size() || _n_speeds != _eff_function_names.size())
102 0 : mooseError("The number of entries of speeds needs to equal the number of entries of "
103 : "Rp_functions and eff_functions");
104 :
105 : // Store functions and check to make sure there is no self-reference.
106 540 : for (unsigned int i = 0; i < _n_speeds; i++)
107 : {
108 478 : if (_Rp_function_names[i] == name() || _eff_function_names[i] == name())
109 0 : mooseError(name(), ": This function cannot use its own name in the 'functions' parameter.");
110 :
111 478 : _Rp_functions[i] = &getFunctionByName(_Rp_function_names[i]);
112 478 : _eff_functions[i] = &getFunctionByName(_eff_function_names[i]);
113 : }
114 62 : }
115 :
116 : void
117 52 : ADShaftConnectedCompressor1PhaseUserObject::initialSetup()
118 : {
119 52 : ADVolumeJunction1PhaseUserObject::initialSetup();
120 :
121 52 : ADShaftConnectableUserObjectInterface::setupConnections(
122 52 : ADVolumeJunctionBaseUserObject::_n_connections, ADVolumeJunctionBaseUserObject::_n_flux_eq);
123 52 : }
124 :
125 : void
126 1639 : ADShaftConnectedCompressor1PhaseUserObject::initialize()
127 : {
128 1639 : ADVolumeJunction1PhaseUserObject::initialize();
129 1639 : ADShaftConnectableUserObjectInterface::initialize();
130 :
131 1639 : _isentropic_torque = 0;
132 1639 : _dissipation_torque = 0;
133 1639 : _friction_torque = 0;
134 1639 : _delta_p = 0;
135 1639 : _Rp = 0;
136 1639 : _eff = 0;
137 1639 : _flow_rel_corr = 0;
138 1639 : _speed_rel_corr = 0;
139 1639 : }
140 :
141 : void
142 2470 : ADShaftConnectedCompressor1PhaseUserObject::execute()
143 : {
144 2470 : ADVolumeJunctionBaseUserObject::storeConnectionData();
145 2470 : ADShaftConnectableUserObjectInterface::setConnectionData(
146 2470 : ADVolumeJunctionBaseUserObject::_flow_channel_dofs);
147 :
148 2470 : const unsigned int c = getBoundaryIDIndex();
149 :
150 2470 : computeFluxesAndResiduals(c);
151 2470 : }
152 :
153 : void
154 2470 : ADShaftConnectedCompressor1PhaseUserObject::computeFluxesAndResiduals(const unsigned int & c)
155 : {
156 2470 : ADVolumeJunction1PhaseUserObject::computeFluxesAndResiduals(c);
157 :
158 : using std::isinf, std::min, std::max, std::abs, std::isnan;
159 :
160 : // inlet c=0 established in component
161 2470 : if (c == 0)
162 : {
163 1235 : const ADReal rhoA_in = _rhoA[0];
164 1235 : const ADReal rhouA_in = _rhouA[0];
165 1235 : const ADReal rhoEA_in = _rhoEA[0];
166 1235 : const ADReal e_in = THM::e_from_arhoA_arhouA_arhoEA(rhoA_in, rhouA_in, rhoEA_in);
167 1235 : const ADReal v_in = _A[0] / rhoA_in;
168 1235 : const ADReal p_in = _fp.p_from_v_e(v_in, e_in);
169 : const ADReal vel_in = rhouA_in / rhoA_in;
170 :
171 : // static entropy is equal to stagnation entropy by definition of the stagnation state
172 1235 : const ADReal s_in = _fp.s_from_v_e(v_in, e_in);
173 1235 : const ADReal s_out = s_in;
174 :
175 : // isentropic: dH/m = Vdp/m
176 : // h0, T0, and c0 are constant for adiabatic flows
177 1235 : const ADReal h0_in = e_in + p_in * v_in + 0.5 * vel_in * vel_in;
178 1235 : const ADReal p0_in = _fp.p_from_h_s(h0_in, s_in);
179 1235 : const ADReal rho0_in = _fp.rho_from_p_s(p0_in, s_in);
180 1235 : const ADReal v0_in = 1.0 / rho0_in;
181 1235 : const ADReal e0_in = _fp.e_from_p_rho(p0_in, rho0_in);
182 1235 : const ADReal c0_in = _fp.c_from_v_e(v0_in, e0_in);
183 1235 : const ADReal flow_A = rhouA_in * _rho0_rated * _c0_rated;
184 :
185 2470 : const ADReal flow_B = _mdot_rated * rho0_in * c0_in;
186 1235 : _flow_rel_corr = flow_A / flow_B;
187 3705 : _speed_rel_corr = (_omega[0] / _omega_rated) * (_c0_rated / c0_in);
188 :
189 : // If _speed_rel_corr == _speeds[0], then the lower index x1 is determined to be -1
190 1235 : if (_speed_rel_corr == _speeds[0])
191 : {
192 41 : _Rp = _Rp_functions[0]->value(_flow_rel_corr, ADPoint());
193 41 : _eff = _eff_functions[0]->value(_flow_rel_corr, ADPoint());
194 : }
195 1194 : else if (isnan(_speed_rel_corr)) // NaN; unguarded, gives segmentation fault
196 : {
197 0 : _Rp = std::nan("");
198 0 : _eff = std::nan("");
199 : }
200 : else // linear interpolation/extrapolation
201 : {
202 : unsigned int x1, x2;
203 1194 : if (_speed_rel_corr < _speeds[0]) // extrapolation past minimum
204 : {
205 : x1 = 0;
206 : x2 = 1;
207 : }
208 939 : else if (_speed_rel_corr > _speeds.back()) // extrapolation past maximum
209 : {
210 0 : x1 = _n_speeds - 1;
211 0 : x2 = x1 - 1;
212 : }
213 : else // interpolation
214 : {
215 939 : auto x1_iter = std::lower_bound(_speeds.begin(), _speeds.end(), _speed_rel_corr);
216 939 : auto x2_iter = std::upper_bound(_speeds.begin(), _speeds.end(), _speed_rel_corr);
217 :
218 939 : x1 = (x1_iter - _speeds.begin()) - 1; // _speeds index for entry <= _speed_rel_corr
219 939 : x2 = (x2_iter - _speeds.begin()); // _speeds index for entry > _speed_rel_corr
220 : }
221 :
222 1194 : const Real x1_spd = _speeds[x1];
223 1194 : const Real x2_spd = _speeds[x2];
224 :
225 1194 : const ADReal y1_Rp = _Rp_functions[x1]->value(_flow_rel_corr, ADPoint());
226 1194 : const ADReal y2_Rp = _Rp_functions[x2]->value(_flow_rel_corr, ADPoint());
227 2388 : const ADReal Rp_m = (y2_Rp - y1_Rp) / (x2_spd - x1_spd);
228 1194 : _Rp = y1_Rp + (_speed_rel_corr - x1_spd) * Rp_m;
229 :
230 1194 : const ADReal y1_eff = _eff_functions[x1]->value(_flow_rel_corr, ADPoint());
231 1194 : const ADReal y2_eff = _eff_functions[x2]->value(_flow_rel_corr, ADPoint());
232 2388 : const ADReal eff_m = (y2_eff - y1_eff) / (x2_spd - x1_spd);
233 1194 : _eff = y1_eff + (_speed_rel_corr - x1_spd) * eff_m;
234 : }
235 :
236 : // Apply bounds
237 1235 : _Rp = max(_Rp_min, min(_Rp_max, _Rp));
238 :
239 : // Invert if treating as turbine
240 : ADReal Rp_comp, eff_comp;
241 1235 : if (_treat_as_turbine)
242 : {
243 510 : Rp_comp = 1.0 / _Rp;
244 510 : eff_comp = 1.0 / _eff;
245 : }
246 : else
247 : {
248 980 : Rp_comp = _Rp;
249 980 : eff_comp = _eff;
250 : }
251 :
252 : const ADReal p0_out = p0_in * Rp_comp;
253 1235 : const ADReal rho0_out_isen = _fp.rho_from_p_s(p0_out, s_out);
254 :
255 1235 : const ADReal e0_out_isen = _fp.e_from_p_rho(p0_out, rho0_out_isen);
256 :
257 1235 : _delta_p = p0_in * (Rp_comp - 1.0);
258 :
259 1235 : if (MooseUtils::absoluteFuzzyEqual(_omega[0], 0.0))
260 : {
261 65 : _isentropic_torque = 0.0;
262 65 : _dissipation_torque = 0.0;
263 : }
264 : else
265 : {
266 1170 : const ADReal h0_out_isen = THM::h_from_e_p_rho(e0_out_isen, p0_out, rho0_out_isen);
267 3510 : _isentropic_torque = -(rhouA_in / _omega[0]) * (h0_out_isen - h0_in); // tau_isen
268 :
269 1170 : const ADReal g_x = h0_out_isen - h0_in + h0_in * eff_comp;
270 : const ADReal h0_out = g_x / eff_comp;
271 :
272 3510 : _dissipation_torque = -(rhouA_in / _omega[0]) * (h0_out - h0_out_isen);
273 : }
274 :
275 1235 : const ADReal alpha = _omega[0] / _omega_rated;
276 :
277 1235 : if (alpha < _speed_cr_I)
278 : {
279 725 : _moment_of_inertia += _inertia_const;
280 : }
281 : else
282 : {
283 1020 : _moment_of_inertia += _inertia_coeff[0] + _inertia_coeff[1] * abs(alpha) +
284 510 : _inertia_coeff[2] * alpha * alpha +
285 1020 : _inertia_coeff[3] * abs(alpha * alpha * alpha);
286 : }
287 :
288 : // friction torque
289 : Real sign;
290 1235 : if (_omega[0] >= 0)
291 1235 : sign = -1;
292 : else
293 0 : sign = 1;
294 1235 : if (alpha < _speed_cr_fr)
295 : {
296 0 : _friction_torque = sign * _tau_fr_const;
297 : }
298 : else
299 : {
300 : _friction_torque =
301 2470 : sign * (_tau_fr_coeff[0] + _tau_fr_coeff[1] * abs(alpha) +
302 3705 : _tau_fr_coeff[2] * alpha * alpha + _tau_fr_coeff[3] * abs(alpha * alpha * alpha));
303 : }
304 :
305 2470 : _torque += _isentropic_torque + _dissipation_torque + _friction_torque;
306 :
307 : // compute momentum and energy source terms
308 : // a negative torque value results in a positive S_energy
309 2470 : const ADReal S_energy = -(_isentropic_torque + _dissipation_torque) * _omega[0];
310 :
311 2470 : const ADRealVectorValue S_momentum = _delta_p * _A_ref * _di_out;
312 :
313 1235 : _residual[VolumeJunction1Phase::RHOEV_INDEX] -= S_energy;
314 :
315 1235 : _residual[VolumeJunction1Phase::RHOUV_INDEX] -= S_momentum(0);
316 1235 : _residual[VolumeJunction1Phase::RHOVV_INDEX] -= S_momentum(1);
317 1235 : _residual[VolumeJunction1Phase::RHOWV_INDEX] -= S_momentum(2);
318 : }
319 2470 : }
320 :
321 : ADReal
322 858 : ADShaftConnectedCompressor1PhaseUserObject::getIsentropicTorque() const
323 : {
324 858 : return _isentropic_torque;
325 : }
326 :
327 : ADReal
328 858 : ADShaftConnectedCompressor1PhaseUserObject::getDissipationTorque() const
329 : {
330 858 : return _dissipation_torque;
331 : }
332 :
333 : ADReal
334 858 : ADShaftConnectedCompressor1PhaseUserObject::getFrictionTorque() const
335 : {
336 858 : return _friction_torque;
337 : }
338 :
339 : ADReal
340 858 : ADShaftConnectedCompressor1PhaseUserObject::getCompressorDeltaP() const
341 : {
342 858 : return _delta_p;
343 : }
344 :
345 : ADReal
346 229 : ADShaftConnectedCompressor1PhaseUserObject::getPressureRatio() const
347 : {
348 229 : return _Rp;
349 : }
350 :
351 : ADReal
352 229 : ADShaftConnectedCompressor1PhaseUserObject::getEfficiency() const
353 : {
354 229 : return _eff;
355 : }
356 : ADReal
357 229 : ADShaftConnectedCompressor1PhaseUserObject::getRelativeCorrectedMassFlowRate() const
358 : {
359 229 : return _flow_rel_corr;
360 : }
361 : ADReal
362 229 : ADShaftConnectedCompressor1PhaseUserObject::getRelativeCorrectedSpeed() const
363 : {
364 229 : return _speed_rel_corr;
365 : }
366 :
367 : void
368 1459 : ADShaftConnectedCompressor1PhaseUserObject::finalize()
369 : {
370 1459 : ADVolumeJunction1PhaseUserObject::finalize();
371 1459 : ADShaftConnectableUserObjectInterface::finalize();
372 :
373 1459 : ADShaftConnectableUserObjectInterface::setupJunctionData(
374 1459 : ADVolumeJunctionBaseUserObject::_scalar_dofs);
375 2918 : ADShaftConnectableUserObjectInterface::setOmegaDofs(getScalarVar("omega", 0));
376 :
377 1459 : comm().sum(_isentropic_torque);
378 1459 : comm().sum(_dissipation_torque);
379 1459 : comm().sum(_friction_torque);
380 1459 : comm().sum(_delta_p);
381 1459 : comm().sum(_Rp);
382 1459 : comm().sum(_eff);
383 1459 : comm().sum(_flow_rel_corr);
384 1459 : comm().sum(_speed_rel_corr);
385 1459 : }
386 :
387 : void
388 180 : ADShaftConnectedCompressor1PhaseUserObject::threadJoin(const UserObject & uo)
389 : {
390 180 : ADVolumeJunction1PhaseUserObject::threadJoin(uo);
391 180 : ADShaftConnectableUserObjectInterface::threadJoin(uo);
392 :
393 : const auto & scpuo = static_cast<const ADShaftConnectedCompressor1PhaseUserObject &>(uo);
394 180 : _isentropic_torque += scpuo._isentropic_torque;
395 180 : _dissipation_torque += scpuo._dissipation_torque;
396 180 : _friction_torque += scpuo._friction_torque;
397 180 : _delta_p += scpuo._delta_p;
398 180 : _Rp += scpuo._Rp;
399 180 : _eff += scpuo._eff;
400 180 : _flow_rel_corr += scpuo._flow_rel_corr;
401 180 : _speed_rel_corr += scpuo._speed_rel_corr;
402 180 : }
|
