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 "ADShaftConnectedTurbine1PhaseUserObject.h"
      11             : #include "ShaftConnectedTurbine1Phase.h"
      12             : #include "SinglePhaseFluidProperties.h"
      13             : #include "VolumeJunction1Phase.h"
      14             : #include "MooseVariableScalar.h"
      15             : #include "THMIndicesVACE.h"
      16             : #include "Function.h"
      17             : #include "Numerics.h"
      18             : #include "metaphysicl/parallel_numberarray.h"
      19             : #include "metaphysicl/parallel_dualnumber.h"
      20             : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
      21             : #include "libmesh/parallel_algebra.h"
      22             : #include "libmesh/utility.h"
      23             : 
      24             : registerMooseObject("ThermalHydraulicsApp", ADShaftConnectedTurbine1PhaseUserObject);
      25             : 
      26             : InputParameters
      27          42 : ADShaftConnectedTurbine1PhaseUserObject::validParams()
      28             : {
      29          42 :   InputParameters params = ADVolumeJunction1PhaseUserObject::validParams();
      30          42 :   params += ADShaftConnectableUserObjectInterface::validParams();
      31             : 
      32          84 :   params.addParam<BoundaryName>("inlet", "Turbine inlet");
      33          84 :   params.addParam<BoundaryName>("outlet", "Turbine outlet");
      34          84 :   params.addRequiredParam<Point>("di_out", "Direction of connected outlet");
      35          84 :   params.addRequiredParam<Real>("omega_rated", "Rated turbine speed [rad/s]");
      36          84 :   params.addRequiredParam<Real>("D_wheel", "Diameter of turbine wheel [m]");
      37          84 :   params.addRequiredParam<Real>("speed_cr_fr", "Turbine speed threshold for friction [-]");
      38          84 :   params.addRequiredParam<Real>("tau_fr_const", "Turbine friction constant [N-m]");
      39          84 :   params.addRequiredParam<std::vector<Real>>("tau_fr_coeff", "Friction coefficients [N-m]");
      40          84 :   params.addRequiredParam<Real>("speed_cr_I", "Turbine speed threshold for inertia [-]");
      41          84 :   params.addRequiredParam<Real>("inertia_const", "Turbine inertia constant [kg-m^2]");
      42          84 :   params.addRequiredParam<std::vector<Real>>("inertia_coeff",
      43             :                                              "Turbine inertia coefficients [kg-m^2]");
      44          84 :   params.addRequiredParam<FunctionName>("head_coefficient",
      45             :                                         "Function to compute data for turbine head [-]");
      46          84 :   params.addRequiredParam<FunctionName>("power_coefficient",
      47             :                                         "Function to compute data for turbine power [-]");
      48          84 :   params.addRequiredParam<std::string>("turbine_name",
      49             :                                        "Name of the instance of this turbine component");
      50          84 :   params.addRequiredCoupledVar("omega", "Shaft rotational speed [rad/s]");
      51             : 
      52          42 :   params.addClassDescription("Computes and caches flux and residual vectors for a 1-phase "
      53             :                              "turbine. Also computes turbine torque "
      54             :                              "and delta_p which is passed to the connected shaft");
      55             : 
      56          42 :   return params;
      57           0 : }
      58             : 
      59          22 : ADShaftConnectedTurbine1PhaseUserObject::ADShaftConnectedTurbine1PhaseUserObject(
      60          22 :     const InputParameters & params)
      61             :   : ADVolumeJunction1PhaseUserObject(params),
      62             :     ADShaftConnectableUserObjectInterface(this),
      63             : 
      64          22 :     _di_out(getParam<Point>("di_out")),
      65          44 :     _omega_rated(getParam<Real>("omega_rated")),
      66          44 :     _D_wheel(getParam<Real>("D_wheel")),
      67          44 :     _speed_cr_fr(getParam<Real>("speed_cr_fr")),
      68          44 :     _tau_fr_const(getParam<Real>("tau_fr_const")),
      69          44 :     _tau_fr_coeff(getParam<std::vector<Real>>("tau_fr_coeff")),
      70          44 :     _speed_cr_I(getParam<Real>("speed_cr_I")),
      71          44 :     _inertia_const(getParam<Real>("inertia_const")),
      72          44 :     _inertia_coeff(getParam<std::vector<Real>>("inertia_coeff")),
      73          22 :     _head_coefficient(getFunction("head_coefficient")),
      74          22 :     _power_coefficient(getFunction("power_coefficient")),
      75          44 :     _turbine_name(getParam<std::string>("turbine_name")),
      76             : 
      77          44 :     _omega(adCoupledScalarValue("omega"))
      78             : {
      79          22 : }
      80             : 
      81             : void
      82          22 : ADShaftConnectedTurbine1PhaseUserObject::initialSetup()
      83             : {
      84          22 :   ADVolumeJunction1PhaseUserObject::initialSetup();
      85             : 
      86          22 :   ADShaftConnectableUserObjectInterface::setupConnections(
      87          22 :       ADVolumeJunctionBaseUserObject::_n_connections, ADVolumeJunctionBaseUserObject::_n_flux_eq);
      88          22 : }
      89             : 
      90             : void
      91        1158 : ADShaftConnectedTurbine1PhaseUserObject::initialize()
      92             : {
      93        1158 :   ADVolumeJunction1PhaseUserObject::initialize();
      94        1158 :   ADShaftConnectableUserObjectInterface::initialize();
      95             : 
      96        1158 :   _driving_torque = 0;
      97        1158 :   _friction_torque = 0;
      98        1158 :   _flow_coeff = 0;
      99        1158 :   _delta_p = 0;
     100        1158 :   _power = 0;
     101        1158 : }
     102             : 
     103             : void
     104        1596 : ADShaftConnectedTurbine1PhaseUserObject::execute()
     105             : {
     106        1596 :   ADVolumeJunctionBaseUserObject::storeConnectionData();
     107        1596 :   ADShaftConnectableUserObjectInterface::setConnectionData(
     108        1596 :       ADVolumeJunctionBaseUserObject::_flow_channel_dofs);
     109             : 
     110        1596 :   const unsigned int c = getBoundaryIDIndex();
     111             : 
     112        1596 :   computeFluxesAndResiduals(c);
     113        1596 : }
     114             : 
     115             : void
     116        1596 : ADShaftConnectedTurbine1PhaseUserObject::computeFluxesAndResiduals(const unsigned int & c)
     117             : {
     118        1596 :   ADVolumeJunction1PhaseUserObject::computeFluxesAndResiduals(c);
     119             : 
     120             :   using std::abs;
     121             : 
     122             :   // inlet c=0 established in component
     123        1596 :   if (c == 0)
     124             :   {
     125             :     const auto & rhoV = _cached_junction_var_values[VolumeJunction1Phase::RHOV_INDEX];
     126             : 
     127         798 :     const ADReal Q_in = (_rhouA[0] / _rhoA[0]) * _A[0];
     128             : 
     129        1596 :     _flow_coeff = Q_in / (_omega[0] * _D_wheel * _D_wheel * _D_wheel);
     130             : 
     131         798 :     const auto head_coeff = _head_coefficient.value(_flow_coeff, ADPoint());
     132             : 
     133         798 :     const ADReal gH = head_coeff * _D_wheel * _D_wheel * _omega[0] * _omega[0];
     134             : 
     135         798 :     const auto power_coeff = _power_coefficient.value(_flow_coeff, ADPoint());
     136             : 
     137             :     // friction torque
     138         798 :     ADReal alpha = _omega[0] / _omega_rated;
     139             :     Real sign;
     140         798 :     if (_omega[0] >= 0)
     141         798 :       sign = -1;
     142             :     else
     143           0 :       sign = 1;
     144         798 :     if (alpha < _speed_cr_fr)
     145             :     {
     146           0 :       _friction_torque = sign * _tau_fr_const;
     147             :     }
     148             :     else
     149             :     {
     150             :       _friction_torque =
     151        1596 :           sign * (_tau_fr_coeff[0] + _tau_fr_coeff[1] * abs(alpha) +
     152        2394 :                   _tau_fr_coeff[2] * alpha * alpha + _tau_fr_coeff[3] * abs(alpha * alpha * alpha));
     153             :     }
     154             : 
     155             :     _driving_torque =
     156        1596 :         power_coeff * (rhoV / _volume) * _omega[0] * _omega[0] * Utility::pow<5>(_D_wheel);
     157             : 
     158        1596 :     _torque += _driving_torque + _friction_torque;
     159             : 
     160         798 :     if (alpha < _speed_cr_I)
     161             :     {
     162         798 :       _moment_of_inertia += _inertia_const;
     163             :     }
     164             :     else
     165             :     {
     166           0 :       _moment_of_inertia += _inertia_coeff[0] + _inertia_coeff[1] * abs(alpha) +
     167           0 :                             _inertia_coeff[2] * alpha * alpha +
     168           0 :                             _inertia_coeff[3] * abs(alpha * alpha * alpha);
     169             :     }
     170             : 
     171             :     // compute momentum and energy source terms
     172             :     // a positive torque value results in a negative S_energy
     173        1596 :     _power = _torque * _omega[0];
     174         798 :     const ADReal S_energy = -_power;
     175             : 
     176        1596 :     _delta_p = (rhoV / _volume) * gH;
     177             : 
     178        1596 :     const ADRealVectorValue S_momentum = -_delta_p * _A_ref * _di_out;
     179             : 
     180         798 :     _residual[VolumeJunction1Phase::RHOEV_INDEX] -= S_energy;
     181             : 
     182         798 :     _residual[VolumeJunction1Phase::RHOUV_INDEX] -= S_momentum(0);
     183         798 :     _residual[VolumeJunction1Phase::RHOVV_INDEX] -= S_momentum(1);
     184         798 :     _residual[VolumeJunction1Phase::RHOWV_INDEX] -= S_momentum(2);
     185             :   }
     186        1596 : }
     187             : 
     188             : ADReal
     189         546 : ADShaftConnectedTurbine1PhaseUserObject::getDrivingTorque() const
     190             : {
     191         546 :   return _driving_torque;
     192             : }
     193             : 
     194             : ADReal
     195         546 : ADShaftConnectedTurbine1PhaseUserObject::getFrictionTorque() const
     196             : {
     197         546 :   return _friction_torque;
     198             : }
     199             : 
     200             : ADReal
     201         546 : ADShaftConnectedTurbine1PhaseUserObject::getFlowCoefficient() const
     202             : {
     203         546 :   return _flow_coeff;
     204             : }
     205             : 
     206             : ADReal
     207         546 : ADShaftConnectedTurbine1PhaseUserObject::getTurbineDeltaP() const
     208             : {
     209         546 :   return _delta_p;
     210             : }
     211             : 
     212             : ADReal
     213         546 : ADShaftConnectedTurbine1PhaseUserObject::getTurbinePower() const
     214             : {
     215         546 :   return _power;
     216             : }
     217             : 
     218             : void
     219        1038 : ADShaftConnectedTurbine1PhaseUserObject::finalize()
     220             : {
     221        1038 :   ADVolumeJunction1PhaseUserObject::finalize();
     222        1038 :   ADShaftConnectableUserObjectInterface::finalize();
     223             : 
     224        1038 :   ADShaftConnectableUserObjectInterface::setupJunctionData(
     225        1038 :       ADVolumeJunctionBaseUserObject::_scalar_dofs);
     226        2076 :   ADShaftConnectableUserObjectInterface::setOmegaDofs(getScalarVar("omega", 0));
     227             : 
     228        1038 :   comm().sum(_driving_torque);
     229        1038 :   comm().sum(_friction_torque);
     230        1038 :   comm().sum(_flow_coeff);
     231        1038 :   comm().sum(_delta_p);
     232        1038 :   comm().sum(_power);
     233        1038 : }
     234             : 
     235             : void
     236         120 : ADShaftConnectedTurbine1PhaseUserObject::threadJoin(const UserObject & uo)
     237             : {
     238         120 :   ADVolumeJunction1PhaseUserObject::threadJoin(uo);
     239         120 :   ADShaftConnectableUserObjectInterface::threadJoin(uo);
     240             : 
     241             :   const auto & scpuo = static_cast<const ADShaftConnectedTurbine1PhaseUserObject &>(uo);
     242         120 :   _driving_torque += scpuo._driving_torque;
     243         120 :   _friction_torque += scpuo._friction_torque;
     244         120 :   _flow_coeff += scpuo._flow_coeff;
     245         120 :   _delta_p += scpuo._delta_p;
     246         120 :   _power += scpuo._power;
     247         120 : }