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 "ADSpecificImpulse1Phase.h"
11 : #include "SinglePhaseFluidProperties.h"
12 : #include "Numerics.h"
13 : #include "THMIndicesVACE.h"
14 : #include "ADBoundaryFluxBase.h"
15 :
16 : registerMooseObject("ThermalHydraulicsApp", ADSpecificImpulse1Phase);
17 :
18 : InputParameters
19 114 : ADSpecificImpulse1Phase::validParams()
20 : {
21 114 : InputParameters params = SidePostprocessor::validParams();
22 228 : params.addRequiredParam<Real>("p_exit", "Outlet pressure at nozzle exit");
23 228 : params.addRequiredParam<UserObjectName>("fp", "Single-phase fluid properties");
24 228 : params.addParam<Real>(
25 228 : "bisection_tolerance", 1e-4, "Tolerance for bisection to find h(entropy_in, p_out)");
26 228 : params.addParam<unsigned int>(
27 : "bisection_max_it",
28 228 : 100,
29 : "Maximum number of iterations for bisection to find h(entropy_in, p_out)");
30 228 : params.addParam<bool>("cumulative",
31 228 : true,
32 : "If specific impulse is accumulated over timesteps. If false, then "
33 : "instantaneous value is computed");
34 228 : params.addCoupledVar("variables", "Single-phase flow variables");
35 684 : params.set<std::vector<VariableName>>("variables") = {"rhoA", "rhouA", "rhoEA", "A"};
36 114 : params.addClassDescription("Estimates specific impulse from fluid state at a boundary");
37 114 : return params;
38 0 : }
39 :
40 42 : ADSpecificImpulse1Phase::ADSpecificImpulse1Phase(const InputParameters & parameters)
41 : : SidePostprocessor(parameters),
42 42 : _n_components(THMVACE1D::N_FLUX_INPUTS),
43 84 : _boundary_name(getParam<std::vector<BoundaryName>>("boundary")[0]),
44 42 : _boundary_uo_name(_boundary_name + ":boundary_uo"),
45 42 : _boundary_uo(getUserObjectByName<ADBoundaryFluxBase>(_boundary_uo_name)),
46 84 : _p_exit(getParam<Real>("p_exit")),
47 84 : _H(getADMaterialPropertyByName<Real>("H")),
48 42 : _v(getADMaterialPropertyByName<Real>("v")),
49 42 : _e(getADMaterialPropertyByName<Real>("e")),
50 42 : _T(getADMaterialPropertyByName<Real>("T")),
51 42 : _fp(getUserObject<SinglePhaseFluidProperties>("fp")),
52 84 : _tol(getParam<Real>("bisection_tolerance")),
53 84 : _max_nit(getParam<unsigned int>("bisection_max_it")),
54 84 : _cumulative(getParam<bool>("cumulative")),
55 84 : _accumulated_mass_flow_rate(declareRestartableData<Real>("accumulated_mass_flow_rate", 0)),
56 84 : _accumulated_thrust(declareRestartableData<Real>("accumulated_thrust", 0)),
57 42 : _value(0.0)
58 : {
59 42 : if (_cumulative && !_fe_problem.isTransient())
60 0 : paramError("cumulative", "Must be false unless problem is transient");
61 :
62 210 : for (unsigned int i = 0; i < _n_components; i++)
63 336 : _U.push_back(&adCoupledValue("variables", i));
64 42 : }
65 :
66 : void
67 150 : ADSpecificImpulse1Phase::threadJoin(const UserObject & y)
68 : {
69 : const ADSpecificImpulse1Phase & pps = static_cast<const ADSpecificImpulse1Phase &>(y);
70 150 : _thrust += pps._thrust;
71 150 : _mass_flow_rate += pps._mass_flow_rate;
72 150 : }
73 :
74 : void
75 800 : ADSpecificImpulse1Phase::initialize()
76 : {
77 800 : _mass_flow_rate = 0;
78 800 : _thrust = 0;
79 800 : }
80 :
81 : Real
82 650 : ADSpecificImpulse1Phase::getValue() const
83 : {
84 650 : return _value;
85 : }
86 :
87 : void
88 650 : ADSpecificImpulse1Phase::finalize()
89 : {
90 650 : gatherSum(_thrust);
91 650 : gatherSum(_mass_flow_rate);
92 :
93 650 : if (_cumulative)
94 : {
95 325 : _accumulated_thrust += _dt * _thrust;
96 325 : _accumulated_mass_flow_rate += _dt * _mass_flow_rate;
97 325 : _value = _accumulated_thrust / _accumulated_mass_flow_rate / THM::gravity_const;
98 : }
99 : else
100 325 : _value = _thrust / _mass_flow_rate / THM::gravity_const;
101 650 : }
102 :
103 : void
104 400 : ADSpecificImpulse1Phase::execute()
105 : {
106 800 : for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
107 : {
108 400 : std::vector<ADReal> U(_n_components);
109 2000 : for (unsigned int i = 0; i < _n_components; i++)
110 1600 : U[i] = (*_U[i])[qp];
111 :
112 400 : const auto & flux = _boundary_uo.getFlux(_current_side, _current_elem->id(), U, _normals[qp]);
113 :
114 400 : _mass_flow_rate += std::abs(MetaPhysicL::raw_value(flux[0]));
115 :
116 : // get entropy at inlet (= entropy at outlet because process is isentropic)
117 : Real entropy_in =
118 400 : _fp.s_from_v_e(MetaPhysicL::raw_value(_v[qp]), MetaPhysicL::raw_value(_e[qp]));
119 :
120 : // compute outlet enthalpy from entropy at inlet (isentropic flow)
121 : // and pressure at outlet, need to do bisection because h_from_s_p does not
122 : // exist, it is better to work with temperature because s(p, T) is available
123 :
124 : // upper bound of temperature is the inlet temperature
125 400 : Real T_up = MetaPhysicL::raw_value(_T[qp]);
126 :
127 : // compute entropy associated with T_up
128 400 : Real entropy_up = _fp.s_from_p_T(_p_exit, T_up);
129 :
130 : // lower bound for temperature is 0.1 K
131 : Real T_low = 0.1;
132 :
133 : // compute entropy associated with T_low
134 400 : Real entropy_low = _fp.s_from_p_T(_p_exit, T_low);
135 :
136 : // compute the midpoint temperature
137 400 : Real T_mid = 0.5 * (T_up + T_low);
138 :
139 : // the current guess for entropy
140 400 : Real entropy_mid = _fp.s_from_p_T(_p_exit, T_mid);
141 :
142 : // main bisection loop
143 : unsigned int nit = 0;
144 4400 : while (std::abs(1 - entropy_mid / entropy_in) > _tol)
145 : {
146 : // check if we are over/underestimating to select either up/down
147 4000 : if ((entropy_mid - entropy_in > 0) == (entropy_up - entropy_low > 0))
148 : {
149 : T_up = T_mid;
150 : entropy_up = entropy_mid;
151 : }
152 : else
153 : {
154 : T_low = T_mid;
155 : entropy_low = entropy_mid;
156 : }
157 :
158 : // update guess for T_mid
159 4000 : T_mid = 0.5 * (T_up + T_low);
160 :
161 : // update the guess for entropy
162 4000 : entropy_mid = _fp.s_from_p_T(_p_exit, T_mid);
163 :
164 4000 : ++nit;
165 :
166 4000 : if (nit == _max_nit)
167 : {
168 0 : mooseDoOnce(mooseWarning("Bisection in ADSpecificImpulse1Phase did not converge"));
169 : break;
170 : }
171 : }
172 :
173 : // the enthalpy evaluated at _p_exit and T_mid
174 400 : Real h_exit = _fp.h_from_p_T(_p_exit, T_mid);
175 :
176 : // compute outlet speed
177 400 : Real vel_exit = std::sqrt(2.0 * (MetaPhysicL::raw_value(_H[qp]) - h_exit));
178 400 : _thrust += std::abs(vel_exit * MetaPhysicL::raw_value(flux[0]));
179 : }
180 400 : }
|
