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 233 : ADShaftConnectedCompressor1PhaseUserObject::validParams()
26 : {
27 233 : InputParameters params = ADVolumeJunction1PhaseUserObject::validParams();
28 233 : params += ADShaftConnectableUserObjectInterface::validParams();
29 :
30 466 : params.addParam<BoundaryName>("inlet", "Compressor inlet");
31 466 : params.addParam<BoundaryName>("outlet", "Compressor outlet");
32 466 : params.addRequiredParam<Point>("di_out", "Direction of connected outlet");
33 466 : params.addRequiredParam<bool>("treat_as_turbine", "Treat the compressor as a turbine?");
34 466 : params.addRequiredParam<Real>("omega_rated", "Rated compressor speed [rad/s]");
35 466 : params.addRequiredParam<Real>("mdot_rated", "Rated compressor mass flow rate [kg/s]");
36 466 : params.addRequiredParam<Real>("rho0_rated",
37 : "Rated compressor inlet stagnation fluid density [kg/m^3]");
38 466 : params.addRequiredParam<Real>("c0_rated", "Rated compressor inlet stagnation sound speed [m/s]");
39 466 : params.addRequiredParam<Real>("speed_cr_fr", "Compressor speed threshold for friction [-]");
40 466 : params.addRequiredParam<Real>("tau_fr_const", "Compressor friction constant [N-m]");
41 466 : params.addRequiredParam<std::vector<Real>>("tau_fr_coeff", "Friction coefficients [N-m]");
42 466 : params.addRequiredParam<Real>("speed_cr_I", "Compressor speed threshold for inertia [-]");
43 466 : params.addRequiredParam<Real>("inertia_const", "Compressor inertia constant [kg-m^2]");
44 466 : params.addRequiredParam<std::vector<Real>>("inertia_coeff",
45 : "Compressor inertia coefficients [kg-m^2]");
46 466 : 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 466 : 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 466 : 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 466 : params.addRequiredParam<Real>("min_pressure_ratio", "Minimum pressure ratio");
61 466 : params.addRequiredParam<Real>("max_pressure_ratio", "Maximum pressure ratio");
62 466 : params.addRequiredParam<std::string>("compressor_name",
63 : "Name of the instance of this compressor component");
64 466 : params.addRequiredCoupledVar("omega", "Shaft rotational speed [rad/s]");
65 :
66 233 : 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 233 : return params;
71 0 : }
72 :
73 125 : ADShaftConnectedCompressor1PhaseUserObject::ADShaftConnectedCompressor1PhaseUserObject(
74 125 : const InputParameters & params)
75 : : ADVolumeJunction1PhaseUserObject(params),
76 : ADShaftConnectableUserObjectInterface(this),
77 :
78 125 : _di_out(getParam<Point>("di_out")),
79 250 : _treat_as_turbine(getParam<bool>("treat_as_turbine")),
80 250 : _omega_rated(getParam<Real>("omega_rated")),
81 250 : _mdot_rated(getParam<Real>("mdot_rated")),
82 250 : _rho0_rated(getParam<Real>("rho0_rated")),
83 250 : _c0_rated(getParam<Real>("c0_rated")),
84 250 : _speed_cr_fr(getParam<Real>("speed_cr_fr")),
85 250 : _tau_fr_const(getParam<Real>("tau_fr_const")),
86 250 : _tau_fr_coeff(getParam<std::vector<Real>>("tau_fr_coeff")),
87 250 : _speed_cr_I(getParam<Real>("speed_cr_I")),
88 250 : _inertia_const(getParam<Real>("inertia_const")),
89 250 : _inertia_coeff(getParam<std::vector<Real>>("inertia_coeff")),
90 250 : _speeds(getParam<std::vector<Real>>("speeds")),
91 250 : _Rp_function_names(getParam<std::vector<FunctionName>>("Rp_functions")),
92 250 : _eff_function_names(getParam<std::vector<FunctionName>>("eff_functions")),
93 125 : _n_speeds(_speeds.size()),
94 125 : _Rp_functions(_n_speeds),
95 125 : _eff_functions(_n_speeds),
96 250 : _Rp_min(getParam<Real>("min_pressure_ratio")),
97 250 : _Rp_max(getParam<Real>("max_pressure_ratio")),
98 250 : _compressor_name(getParam<std::string>("compressor_name")),
99 250 : _omega(adCoupledScalarValue("omega"))
100 : {
101 125 : 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 1092 : for (unsigned int i = 0; i < _n_speeds; i++)
107 : {
108 967 : 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 967 : _Rp_functions[i] = &getFunctionByName(_Rp_function_names[i]);
112 967 : _eff_functions[i] = &getFunctionByName(_eff_function_names[i]);
113 : }
114 125 : }
115 :
116 : void
117 109 : ADShaftConnectedCompressor1PhaseUserObject::initialSetup()
118 : {
119 109 : ADVolumeJunction1PhaseUserObject::initialSetup();
120 :
121 109 : ADShaftConnectableUserObjectInterface::setupConnections(
122 109 : ADVolumeJunctionBaseUserObject::_n_connections, ADVolumeJunctionBaseUserObject::_n_flux_eq);
123 109 : }
124 :
125 : void
126 12847 : ADShaftConnectedCompressor1PhaseUserObject::initialize()
127 : {
128 12847 : ADVolumeJunction1PhaseUserObject::initialize();
129 12847 : ADShaftConnectableUserObjectInterface::initialize();
130 :
131 12847 : _isentropic_torque = 0;
132 12847 : _dissipation_torque = 0;
133 12847 : _friction_torque = 0;
134 12847 : _delta_p = 0;
135 12847 : _Rp = 0;
136 12847 : _eff = 0;
137 12847 : _flow_rel_corr = 0;
138 12847 : _speed_rel_corr = 0;
139 12847 : }
140 :
141 : void
142 21358 : ADShaftConnectedCompressor1PhaseUserObject::execute()
143 : {
144 21358 : ADVolumeJunctionBaseUserObject::storeConnectionData();
145 21358 : ADShaftConnectableUserObjectInterface::setConnectionData(
146 21358 : ADVolumeJunctionBaseUserObject::_flow_channel_dofs);
147 :
148 21358 : const unsigned int c = getBoundaryIDIndex();
149 :
150 21358 : computeFluxesAndResiduals(c);
151 21358 : }
152 :
153 : void
154 21358 : ADShaftConnectedCompressor1PhaseUserObject::computeFluxesAndResiduals(const unsigned int & c)
155 : {
156 21358 : 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 21358 : if (c == 0)
162 : {
163 10679 : const ADReal rhoA_in = _rhoA[0];
164 10679 : const ADReal rhouA_in = _rhouA[0];
165 10679 : const ADReal rhoEA_in = _rhoEA[0];
166 10679 : const ADReal e_in = THM::e_from_arhoA_arhouA_arhoEA(rhoA_in, rhouA_in, rhoEA_in);
167 10679 : const ADReal v_in = _A[0] / rhoA_in;
168 10679 : 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 10679 : const ADReal s_in = _fp.s_from_v_e(v_in, e_in);
173 10679 : const ADReal s_out = s_in;
174 :
175 : // isentropic: dH/m = Vdp/m
176 : // h0, T0, and c0 are constant for adiabatic flows
177 10679 : const ADReal h0_in = e_in + p_in * v_in + 0.5 * vel_in * vel_in;
178 10679 : const ADReal p0_in = _fp.p_from_h_s(h0_in, s_in);
179 10679 : const ADReal rho0_in = _fp.rho_from_p_s(p0_in, s_in);
180 10679 : const ADReal v0_in = 1.0 / rho0_in;
181 10679 : const ADReal e0_in = _fp.e_from_p_rho(p0_in, rho0_in);
182 10679 : const ADReal c0_in = _fp.c_from_v_e(v0_in, e0_in);
183 10679 : const ADReal flow_A = rhouA_in * _rho0_rated * _c0_rated;
184 :
185 21358 : const ADReal flow_B = _mdot_rated * rho0_in * c0_in;
186 10679 : _flow_rel_corr = flow_A / flow_B;
187 32037 : _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 10679 : if (_speed_rel_corr == _speeds[0])
191 : {
192 68 : _Rp = _Rp_functions[0]->value(_flow_rel_corr, ADPoint());
193 68 : _eff = _eff_functions[0]->value(_flow_rel_corr, ADPoint());
194 : }
195 10611 : 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 10611 : if (_speed_rel_corr < _speeds[0]) // extrapolation past minimum
204 : {
205 : x1 = 0;
206 : x2 = 1;
207 : }
208 5819 : 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 5819 : auto x1_iter = std::lower_bound(_speeds.begin(), _speeds.end(), _speed_rel_corr);
216 5819 : auto x2_iter = std::upper_bound(_speeds.begin(), _speeds.end(), _speed_rel_corr);
217 :
218 5819 : x1 = (x1_iter - _speeds.begin()) - 1; // _speeds index for entry <= _speed_rel_corr
219 5819 : x2 = (x2_iter - _speeds.begin()); // _speeds index for entry > _speed_rel_corr
220 : }
221 :
222 10611 : const Real x1_spd = _speeds[x1];
223 10611 : const Real x2_spd = _speeds[x2];
224 :
225 10611 : const ADReal y1_Rp = _Rp_functions[x1]->value(_flow_rel_corr, ADPoint());
226 10611 : const ADReal y2_Rp = _Rp_functions[x2]->value(_flow_rel_corr, ADPoint());
227 21222 : const ADReal Rp_m = (y2_Rp - y1_Rp) / (x2_spd - x1_spd);
228 10611 : _Rp = y1_Rp + (_speed_rel_corr - x1_spd) * Rp_m;
229 :
230 10611 : const ADReal y1_eff = _eff_functions[x1]->value(_flow_rel_corr, ADPoint());
231 10611 : const ADReal y2_eff = _eff_functions[x2]->value(_flow_rel_corr, ADPoint());
232 21222 : const ADReal eff_m = (y2_eff - y1_eff) / (x2_spd - x1_spd);
233 10611 : _eff = y1_eff + (_speed_rel_corr - x1_spd) * eff_m;
234 : }
235 :
236 : // Apply bounds
237 10679 : _Rp = max(_Rp_min, min(_Rp_max, _Rp));
238 :
239 : // Invert if treating as turbine
240 : ADReal Rp_comp, eff_comp;
241 10679 : if (_treat_as_turbine)
242 : {
243 9584 : Rp_comp = 1.0 / _Rp;
244 9584 : eff_comp = 1.0 / _eff;
245 : }
246 : else
247 : {
248 5887 : Rp_comp = _Rp;
249 5887 : eff_comp = _eff;
250 : }
251 :
252 : const ADReal p0_out = p0_in * Rp_comp;
253 10679 : const ADReal rho0_out_isen = _fp.rho_from_p_s(p0_out, s_out);
254 :
255 10679 : const ADReal e0_out_isen = _fp.e_from_p_rho(p0_out, rho0_out_isen);
256 :
257 10679 : _delta_p = p0_in * (Rp_comp - 1.0);
258 :
259 10679 : if (MooseUtils::absoluteFuzzyEqual(_omega[0], 0.0))
260 : {
261 116 : _isentropic_torque = 0.0;
262 116 : _dissipation_torque = 0.0;
263 : }
264 : else
265 : {
266 10563 : const ADReal h0_out_isen = THM::h_from_e_p_rho(e0_out_isen, p0_out, rho0_out_isen);
267 31689 : _isentropic_torque = -(rhouA_in / _omega[0]) * (h0_out_isen - h0_in); // tau_isen
268 :
269 10563 : const ADReal g_x = h0_out_isen - h0_in + h0_in * eff_comp;
270 : const ADReal h0_out = g_x / eff_comp;
271 :
272 31689 : _dissipation_torque = -(rhouA_in / _omega[0]) * (h0_out - h0_out_isen);
273 : }
274 :
275 10679 : const ADReal alpha = _omega[0] / _omega_rated;
276 :
277 10679 : if (alpha < _speed_cr_I)
278 : {
279 1095 : _moment_of_inertia += _inertia_const;
280 : }
281 : else
282 : {
283 19168 : _moment_of_inertia += _inertia_coeff[0] + _inertia_coeff[1] * abs(alpha) +
284 9584 : _inertia_coeff[2] * alpha * alpha +
285 19168 : _inertia_coeff[3] * abs(alpha * alpha * alpha);
286 : }
287 :
288 : // friction torque
289 : Real sign;
290 10679 : if (_omega[0] >= 0)
291 10679 : sign = -1;
292 : else
293 0 : sign = 1;
294 10679 : if (alpha < _speed_cr_fr)
295 : {
296 0 : _friction_torque = sign * _tau_fr_const;
297 : }
298 : else
299 : {
300 : _friction_torque =
301 21358 : sign * (_tau_fr_coeff[0] + _tau_fr_coeff[1] * abs(alpha) +
302 32037 : _tau_fr_coeff[2] * alpha * alpha + _tau_fr_coeff[3] * abs(alpha * alpha * alpha));
303 : }
304 :
305 21358 : _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 21358 : const ADReal S_energy = -(_isentropic_torque + _dissipation_torque) * _omega[0];
310 :
311 21358 : const ADRealVectorValue S_momentum = _delta_p * _A_ref * _di_out;
312 :
313 10679 : _residual[VolumeJunction1Phase::RHOEV_INDEX] -= S_energy;
314 :
315 10679 : _residual[VolumeJunction1Phase::RHOUV_INDEX] -= S_momentum(0);
316 10679 : _residual[VolumeJunction1Phase::RHOVV_INDEX] -= S_momentum(1);
317 10679 : _residual[VolumeJunction1Phase::RHOWV_INDEX] -= S_momentum(2);
318 : }
319 21358 : }
320 :
321 : ADReal
322 7082 : ADShaftConnectedCompressor1PhaseUserObject::getIsentropicTorque() const
323 : {
324 7082 : return _isentropic_torque;
325 : }
326 :
327 : ADReal
328 7082 : ADShaftConnectedCompressor1PhaseUserObject::getDissipationTorque() const
329 : {
330 7082 : return _dissipation_torque;
331 : }
332 :
333 : ADReal
334 7082 : ADShaftConnectedCompressor1PhaseUserObject::getFrictionTorque() const
335 : {
336 7082 : return _friction_torque;
337 : }
338 :
339 : ADReal
340 7082 : ADShaftConnectedCompressor1PhaseUserObject::getCompressorDeltaP() const
341 : {
342 7082 : return _delta_p;
343 : }
344 :
345 : ADReal
346 1214 : ADShaftConnectedCompressor1PhaseUserObject::getPressureRatio() const
347 : {
348 1214 : return _Rp;
349 : }
350 :
351 : ADReal
352 1214 : ADShaftConnectedCompressor1PhaseUserObject::getEfficiency() const
353 : {
354 1214 : return _eff;
355 : }
356 : ADReal
357 1214 : ADShaftConnectedCompressor1PhaseUserObject::getRelativeCorrectedMassFlowRate() const
358 : {
359 1214 : return _flow_rel_corr;
360 : }
361 : ADReal
362 1214 : ADShaftConnectedCompressor1PhaseUserObject::getRelativeCorrectedSpeed() const
363 : {
364 1214 : return _speed_rel_corr;
365 : }
366 :
367 : void
368 11219 : ADShaftConnectedCompressor1PhaseUserObject::finalize()
369 : {
370 11219 : ADVolumeJunction1PhaseUserObject::finalize();
371 11219 : ADShaftConnectableUserObjectInterface::finalize();
372 :
373 11219 : ADShaftConnectableUserObjectInterface::setupJunctionData(
374 11219 : ADVolumeJunctionBaseUserObject::_scalar_dofs);
375 22438 : ADShaftConnectableUserObjectInterface::setOmegaDofs(getScalarVar("omega", 0));
376 :
377 11219 : comm().sum(_isentropic_torque);
378 11219 : comm().sum(_dissipation_torque);
379 11219 : comm().sum(_friction_torque);
380 11219 : comm().sum(_delta_p);
381 11219 : comm().sum(_Rp);
382 11219 : comm().sum(_eff);
383 11219 : comm().sum(_flow_rel_corr);
384 11219 : comm().sum(_speed_rel_corr);
385 11219 : }
386 :
387 : void
388 1628 : ADShaftConnectedCompressor1PhaseUserObject::threadJoin(const UserObject & uo)
389 : {
390 1628 : ADVolumeJunction1PhaseUserObject::threadJoin(uo);
391 1628 : ADShaftConnectableUserObjectInterface::threadJoin(uo);
392 :
393 : const auto & scpuo = static_cast<const ADShaftConnectedCompressor1PhaseUserObject &>(uo);
394 1628 : _isentropic_torque += scpuo._isentropic_torque;
395 1628 : _dissipation_torque += scpuo._dissipation_torque;
396 1628 : _friction_torque += scpuo._friction_torque;
397 1628 : _delta_p += scpuo._delta_p;
398 1628 : _Rp += scpuo._Rp;
399 1628 : _eff += scpuo._eff;
400 1628 : _flow_rel_corr += scpuo._flow_rel_corr;
401 1628 : _speed_rel_corr += scpuo._speed_rel_corr;
402 1628 : }
|
