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