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 : // inlet c=0 established in component
159 21358 : if (c == 0)
160 : {
161 10679 : const ADReal rhoA_in = _rhoA[0];
162 10679 : const ADReal rhouA_in = _rhouA[0];
163 10679 : const ADReal rhoEA_in = _rhoEA[0];
164 10679 : const ADReal e_in = THM::e_from_arhoA_arhouA_arhoEA(rhoA_in, rhouA_in, rhoEA_in);
165 10679 : const ADReal v_in = _A[0] / rhoA_in;
166 10679 : const ADReal p_in = _fp.p_from_v_e(v_in, e_in);
167 : const ADReal vel_in = rhouA_in / rhoA_in;
168 :
169 : // static entropy is equal to stagnation entropy by definition of the stagnation state
170 10679 : const ADReal s_in = _fp.s_from_v_e(v_in, e_in);
171 10679 : const ADReal s_out = s_in;
172 :
173 : // isentropic: dH/m = Vdp/m
174 : // h0, T0, and c0 are constant for adiabatic flows
175 10679 : const ADReal h0_in = e_in + p_in * v_in + 0.5 * vel_in * vel_in;
176 10679 : const ADReal p0_in = _fp.p_from_h_s(h0_in, s_in);
177 10679 : const ADReal rho0_in = _fp.rho_from_p_s(p0_in, s_in);
178 10679 : const ADReal v0_in = 1.0 / rho0_in;
179 10679 : const ADReal e0_in = _fp.e_from_p_rho(p0_in, rho0_in);
180 10679 : const ADReal c0_in = _fp.c_from_v_e(v0_in, e0_in);
181 10679 : const ADReal flow_A = rhouA_in * _rho0_rated * _c0_rated;
182 :
183 21358 : const ADReal flow_B = _mdot_rated * rho0_in * c0_in;
184 10679 : _flow_rel_corr = flow_A / flow_B;
185 32037 : _speed_rel_corr = (_omega[0] / _omega_rated) * (_c0_rated / c0_in);
186 :
187 : // If _speed_rel_corr == _speeds[0], then the lower index x1 is determined to be -1
188 10679 : if (_speed_rel_corr == _speeds[0])
189 : {
190 68 : _Rp = _Rp_functions[0]->value(_flow_rel_corr, ADPoint());
191 68 : _eff = _eff_functions[0]->value(_flow_rel_corr, ADPoint());
192 : }
193 10611 : else if (std::isnan(_speed_rel_corr)) // NaN; unguarded, gives segmentation fault
194 : {
195 0 : _Rp = std::nan("");
196 0 : _eff = std::nan("");
197 : }
198 : else // linear interpolation/extrapolation
199 : {
200 : unsigned int x1, x2;
201 10611 : if (_speed_rel_corr < _speeds[0]) // extrapolation past minimum
202 : {
203 : x1 = 0;
204 : x2 = 1;
205 : }
206 5819 : else if (_speed_rel_corr > _speeds.back()) // extrapolation past maximum
207 : {
208 0 : x1 = _n_speeds - 1;
209 0 : x2 = x1 - 1;
210 : }
211 : else // interpolation
212 : {
213 5819 : auto x1_iter = std::lower_bound(_speeds.begin(), _speeds.end(), _speed_rel_corr);
214 5819 : auto x2_iter = std::upper_bound(_speeds.begin(), _speeds.end(), _speed_rel_corr);
215 :
216 5819 : x1 = (x1_iter - _speeds.begin()) - 1; // _speeds index for entry <= _speed_rel_corr
217 5819 : x2 = (x2_iter - _speeds.begin()); // _speeds index for entry > _speed_rel_corr
218 : }
219 :
220 10611 : const Real x1_spd = _speeds[x1];
221 10611 : const Real x2_spd = _speeds[x2];
222 :
223 10611 : const ADReal y1_Rp = _Rp_functions[x1]->value(_flow_rel_corr, ADPoint());
224 10611 : const ADReal y2_Rp = _Rp_functions[x2]->value(_flow_rel_corr, ADPoint());
225 21222 : const ADReal Rp_m = (y2_Rp - y1_Rp) / (x2_spd - x1_spd);
226 10611 : _Rp = y1_Rp + (_speed_rel_corr - x1_spd) * Rp_m;
227 :
228 10611 : const ADReal y1_eff = _eff_functions[x1]->value(_flow_rel_corr, ADPoint());
229 10611 : const ADReal y2_eff = _eff_functions[x2]->value(_flow_rel_corr, ADPoint());
230 21222 : const ADReal eff_m = (y2_eff - y1_eff) / (x2_spd - x1_spd);
231 10611 : _eff = y1_eff + (_speed_rel_corr - x1_spd) * eff_m;
232 : }
233 :
234 : // Apply bounds
235 10679 : _Rp = std::max(_Rp_min, std::min(_Rp_max, _Rp));
236 :
237 : // Invert if treating as turbine
238 : ADReal Rp_comp, eff_comp;
239 10679 : if (_treat_as_turbine)
240 : {
241 9584 : Rp_comp = 1.0 / _Rp;
242 9584 : eff_comp = 1.0 / _eff;
243 : }
244 : else
245 : {
246 5887 : Rp_comp = _Rp;
247 5887 : eff_comp = _eff;
248 : }
249 :
250 : const ADReal p0_out = p0_in * Rp_comp;
251 10679 : const ADReal rho0_out_isen = _fp.rho_from_p_s(p0_out, s_out);
252 :
253 10679 : const ADReal e0_out_isen = _fp.e_from_p_rho(p0_out, rho0_out_isen);
254 :
255 10679 : _delta_p = p0_in * (Rp_comp - 1.0);
256 :
257 10679 : if (MooseUtils::absoluteFuzzyEqual(_omega[0], 0.0))
258 : {
259 116 : _isentropic_torque = 0.0;
260 116 : _dissipation_torque = 0.0;
261 : }
262 : else
263 : {
264 10563 : const ADReal h0_out_isen = THM::h_from_e_p_rho(e0_out_isen, p0_out, rho0_out_isen);
265 31689 : _isentropic_torque = -(rhouA_in / _omega[0]) * (h0_out_isen - h0_in); // tau_isen
266 :
267 10563 : const ADReal g_x = h0_out_isen - h0_in + h0_in * eff_comp;
268 : const ADReal h0_out = g_x / eff_comp;
269 :
270 31689 : _dissipation_torque = -(rhouA_in / _omega[0]) * (h0_out - h0_out_isen);
271 : }
272 :
273 10679 : const ADReal alpha = _omega[0] / _omega_rated;
274 :
275 10679 : if (alpha < _speed_cr_I)
276 : {
277 1095 : _moment_of_inertia += _inertia_const;
278 : }
279 : else
280 : {
281 19168 : _moment_of_inertia += _inertia_coeff[0] + _inertia_coeff[1] * std::abs(alpha) +
282 9584 : _inertia_coeff[2] * alpha * alpha +
283 19168 : _inertia_coeff[3] * std::abs(alpha * alpha * alpha);
284 : }
285 :
286 : // friction torque
287 : Real sign;
288 10679 : if (_omega[0] >= 0)
289 10679 : sign = -1;
290 : else
291 0 : sign = 1;
292 10679 : if (alpha < _speed_cr_fr)
293 : {
294 0 : _friction_torque = sign * _tau_fr_const;
295 : }
296 : else
297 : {
298 21358 : _friction_torque = sign * (_tau_fr_coeff[0] + _tau_fr_coeff[1] * std::abs(alpha) +
299 10679 : _tau_fr_coeff[2] * alpha * alpha +
300 21358 : _tau_fr_coeff[3] * std::abs(alpha * alpha * alpha));
301 : }
302 :
303 21358 : _torque += _isentropic_torque + _dissipation_torque + _friction_torque;
304 :
305 : // compute momentum and energy source terms
306 : // a negative torque value results in a positive S_energy
307 21358 : const ADReal S_energy = -(_isentropic_torque + _dissipation_torque) * _omega[0];
308 :
309 21358 : const ADRealVectorValue S_momentum = _delta_p * _A_ref * _di_out;
310 :
311 10679 : _residual[VolumeJunction1Phase::RHOEV_INDEX] -= S_energy;
312 :
313 10679 : _residual[VolumeJunction1Phase::RHOUV_INDEX] -= S_momentum(0);
314 10679 : _residual[VolumeJunction1Phase::RHOVV_INDEX] -= S_momentum(1);
315 10679 : _residual[VolumeJunction1Phase::RHOWV_INDEX] -= S_momentum(2);
316 : }
317 21358 : }
318 :
319 : ADReal
320 7082 : ADShaftConnectedCompressor1PhaseUserObject::getIsentropicTorque() const
321 : {
322 7082 : return _isentropic_torque;
323 : }
324 :
325 : ADReal
326 7082 : ADShaftConnectedCompressor1PhaseUserObject::getDissipationTorque() const
327 : {
328 7082 : return _dissipation_torque;
329 : }
330 :
331 : ADReal
332 7082 : ADShaftConnectedCompressor1PhaseUserObject::getFrictionTorque() const
333 : {
334 7082 : return _friction_torque;
335 : }
336 :
337 : ADReal
338 7082 : ADShaftConnectedCompressor1PhaseUserObject::getCompressorDeltaP() const
339 : {
340 7082 : return _delta_p;
341 : }
342 :
343 : ADReal
344 1214 : ADShaftConnectedCompressor1PhaseUserObject::getPressureRatio() const
345 : {
346 1214 : return _Rp;
347 : }
348 :
349 : ADReal
350 1214 : ADShaftConnectedCompressor1PhaseUserObject::getEfficiency() const
351 : {
352 1214 : return _eff;
353 : }
354 : ADReal
355 1214 : ADShaftConnectedCompressor1PhaseUserObject::getRelativeCorrectedMassFlowRate() const
356 : {
357 1214 : return _flow_rel_corr;
358 : }
359 : ADReal
360 1214 : ADShaftConnectedCompressor1PhaseUserObject::getRelativeCorrectedSpeed() const
361 : {
362 1214 : return _speed_rel_corr;
363 : }
364 :
365 : void
366 11219 : ADShaftConnectedCompressor1PhaseUserObject::finalize()
367 : {
368 11219 : ADVolumeJunction1PhaseUserObject::finalize();
369 11219 : ADShaftConnectableUserObjectInterface::finalize();
370 :
371 11219 : ADShaftConnectableUserObjectInterface::setupJunctionData(
372 11219 : ADVolumeJunctionBaseUserObject::_scalar_dofs);
373 22438 : ADShaftConnectableUserObjectInterface::setOmegaDofs(getScalarVar("omega", 0));
374 :
375 11219 : comm().sum(_isentropic_torque);
376 11219 : comm().sum(_dissipation_torque);
377 11219 : comm().sum(_friction_torque);
378 11219 : comm().sum(_delta_p);
379 11219 : comm().sum(_Rp);
380 11219 : comm().sum(_eff);
381 11219 : comm().sum(_flow_rel_corr);
382 11219 : comm().sum(_speed_rel_corr);
383 11219 : }
384 :
385 : void
386 1628 : ADShaftConnectedCompressor1PhaseUserObject::threadJoin(const UserObject & uo)
387 : {
388 1628 : ADVolumeJunction1PhaseUserObject::threadJoin(uo);
389 1628 : ADShaftConnectableUserObjectInterface::threadJoin(uo);
390 :
391 : const auto & scpuo = static_cast<const ADShaftConnectedCompressor1PhaseUserObject &>(uo);
392 1628 : _isentropic_torque += scpuo._isentropic_torque;
393 1628 : _dissipation_torque += scpuo._dissipation_torque;
394 1628 : _friction_torque += scpuo._friction_torque;
395 1628 : _delta_p += scpuo._delta_p;
396 1628 : _Rp += scpuo._Rp;
397 1628 : _eff += scpuo._eff;
398 1628 : _flow_rel_corr += scpuo._flow_rel_corr;
399 1628 : _speed_rel_corr += scpuo._speed_rel_corr;
400 1628 : }
|
