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