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 "ADVolumeJunction1PhaseUserObject.h"
      11             : #include "SinglePhaseFluidProperties.h"
      12             : #include "THMIndicesVACE.h"
      13             : #include "VolumeJunction1Phase.h"
      14             : #include "ADNumericalFlux3EqnBase.h"
      15             : #include "Numerics.h"
      16             : #include "THMUtils.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", ADVolumeJunction1PhaseUserObject);
      23             : 
      24             : InputParameters
      25        1205 : ADVolumeJunction1PhaseUserObject::validParams()
      26             : {
      27        1205 :   InputParameters params = ADVolumeJunctionBaseUserObject::validParams();
      28             : 
      29        2410 :   params.addRequiredCoupledVar("A", "Cross-sectional area of connected flow channels");
      30        2410 :   params.addRequiredCoupledVar("rhoA", "rho*A of the connected flow channels");
      31        2410 :   params.addRequiredCoupledVar("rhouA", "rhou*A of the connected flow channels");
      32        2410 :   params.addRequiredCoupledVar("rhoEA", "rhoE*A of the connected flow channels");
      33             : 
      34        2410 :   params.addRequiredCoupledVar("rhoV", "rho*V of the junction");
      35        2410 :   params.addRequiredCoupledVar("rhouV", "rho*u*V of the junction");
      36        2410 :   params.addRequiredCoupledVar("rhovV", "rho*v*V of the junction");
      37        2410 :   params.addRequiredCoupledVar("rhowV", "rho*w*V of the junction");
      38        2410 :   params.addRequiredCoupledVar("rhoEV", "rho*E*V of the junction");
      39             : 
      40        2410 :   params.addRequiredParam<UserObjectName>("fp", "Fluid properties user object name");
      41             : 
      42        2410 :   params.addRequiredParam<Real>("K", "Form loss coefficient [-]");
      43        2410 :   params.addRequiredParam<Real>("A_ref", "Reference area [m^2]");
      44             : 
      45        2410 :   params.addParam<bool>("apply_velocity_scaling",
      46        2410 :                         false,
      47             :                         "Set to true to apply the scaling to the normal velocity. See "
      48             :                         "documentation for more information.");
      49             : 
      50        2410 :   params.addClassDescription(
      51             :       "Computes and caches flux and residual vectors for a 1-phase volume junction");
      52             : 
      53        2410 :   params.declareControllable("K");
      54        1205 :   return params;
      55           0 : }
      56             : 
      57         641 : ADVolumeJunction1PhaseUserObject::ADVolumeJunction1PhaseUserObject(const InputParameters & params)
      58             :   : ADVolumeJunctionBaseUserObject(params),
      59             : 
      60         641 :     _A(adCoupledValue("A")),
      61         641 :     _rhoA(adCoupledValue("rhoA")),
      62         641 :     _rhouA(adCoupledValue("rhouA")),
      63         641 :     _rhoEA(adCoupledValue("rhoEA")),
      64             : 
      65        1282 :     _K(getParam<Real>("K")),
      66        1282 :     _A_ref(getParam<Real>("A_ref")),
      67             : 
      68        1282 :     _apply_velocity_scaling(getParam<bool>("apply_velocity_scaling")),
      69             : 
      70        1282 :     _fp(getUserObject<SinglePhaseFluidProperties>("fp"))
      71             : {
      72         641 :   _flow_variable_names.resize(THMVACE1D::N_FLUX_OUTPUTS);
      73             :   _flow_variable_names[THMVACE1D::RHOA] = "rhoA";
      74             :   _flow_variable_names[THMVACE1D::RHOUA] = "rhouA";
      75             :   _flow_variable_names[THMVACE1D::RHOEA] = "rhoEA";
      76             : 
      77         641 :   _scalar_variable_names.resize(VolumeJunction1Phase::N_EQ);
      78             :   _scalar_variable_names[VolumeJunction1Phase::RHOV_INDEX] = "rhoV";
      79             :   _scalar_variable_names[VolumeJunction1Phase::RHOUV_INDEX] = "rhouV";
      80             :   _scalar_variable_names[VolumeJunction1Phase::RHOVV_INDEX] = "rhovV";
      81             :   _scalar_variable_names[VolumeJunction1Phase::RHOWV_INDEX] = "rhowV";
      82             :   _scalar_variable_names[VolumeJunction1Phase::RHOEV_INDEX] = "rhoEV";
      83             : 
      84         641 :   _junction_var_values.resize(VolumeJunction1Phase::N_EQ);
      85         641 :   _junction_var_values[VolumeJunction1Phase::RHOV_INDEX] = &coupledJunctionValue("rhoV");
      86         641 :   _junction_var_values[VolumeJunction1Phase::RHOUV_INDEX] = &coupledJunctionValue("rhouV");
      87         641 :   _junction_var_values[VolumeJunction1Phase::RHOVV_INDEX] = &coupledJunctionValue("rhovV");
      88         641 :   _junction_var_values[VolumeJunction1Phase::RHOWV_INDEX] = &coupledJunctionValue("rhowV");
      89         641 :   _junction_var_values[VolumeJunction1Phase::RHOEV_INDEX] = &coupledJunctionValue("rhoEV");
      90             : 
      91         641 :   _numerical_flux_uo.resize(_n_connections);
      92        2052 :   for (std::size_t i = 0; i < _n_connections; i++)
      93        1411 :     _numerical_flux_uo[i] = &getUserObjectByName<ADNumericalFlux3EqnBase>(_numerical_flux_names[i]);
      94         641 : }
      95             : 
      96             : void
      97       66856 : ADVolumeJunction1PhaseUserObject::computeFluxesAndResiduals(const unsigned int & c)
      98             : {
      99       66856 :   const Real din = _normal[c];
     100       66856 :   const RealVectorValue di = _dir[0];
     101             :   const RealVectorValue ni = di * din;
     102       66856 :   const Real nJi_dot_di = -din;
     103             : 
     104       66856 :   std::vector<ADReal> Ui(THMVACE3D::N_FLUX_INPUTS, 0.);
     105       66856 :   Ui[THMVACE3D::RHOA] = _rhoA[0];
     106      133712 :   Ui[THMVACE3D::RHOUA] = _rhouA[0] * di(0);
     107       66856 :   Ui[THMVACE3D::RHOVA] = _rhouA[0] * di(1);
     108       66856 :   Ui[THMVACE3D::RHOWA] = _rhouA[0] * di(2);
     109       66856 :   Ui[THMVACE3D::RHOEA] = _rhoEA[0];
     110       66856 :   Ui[THMVACE3D::AREA] = _A[0];
     111             : 
     112       66856 :   const auto flux_3d = compute3DFlux(*(_numerical_flux_uo[c]), Ui, ni);
     113             : 
     114       66856 :   _flux[c].resize(THMVACE1D::N_FLUX_OUTPUTS);
     115       66856 :   _flux[c][THMVACE1D::MASS] = flux_3d[THMVACE3D::MASS] * nJi_dot_di;
     116       66856 :   _flux[c][THMVACE1D::MOMENTUM] = flux_3d[THMVACE3D::MOM_NORM];
     117       66856 :   _flux[c][THMVACE1D::ENERGY] = flux_3d[THMVACE3D::ENERGY] * nJi_dot_di;
     118             : 
     119       66856 :   const bool is_primary_connection = (c == 0);
     120       66856 :   const auto residual = computeResidual(flux_3d, Ui, ni, is_primary_connection);
     121             : 
     122      401136 :   for (const auto i : index_range(_residual))
     123      334280 :     _residual[i] += residual[i];
     124       66856 : }
     125             : 
     126             : std::vector<ADReal>
     127       69679 : ADVolumeJunction1PhaseUserObject::compute3DFlux(const ADNumericalFlux3EqnBase & numerical_flux,
     128             :                                                 const std::vector<ADReal> & Ui,
     129             :                                                 const RealVectorValue & ni) const
     130             : {
     131             :   using std::abs, std::max, std::min;
     132             : 
     133             :   const RealVectorValue nJi = -ni;
     134             : 
     135             :   RealVectorValue t1, t2;
     136       69679 :   THM::computeOrthogonalDirections(nJi, t1, t2);
     137             : 
     138             :   const auto & Ai = Ui[THMVACE3D::AREA];
     139             : 
     140             :   const auto & rhoV = _cached_junction_var_values[VolumeJunction1Phase::RHOV_INDEX];
     141             :   const auto & rhouV = _cached_junction_var_values[VolumeJunction1Phase::RHOUV_INDEX];
     142             :   const auto & rhovV = _cached_junction_var_values[VolumeJunction1Phase::RHOVV_INDEX];
     143             :   const auto & rhowV = _cached_junction_var_values[VolumeJunction1Phase::RHOWV_INDEX];
     144             :   const auto & rhoEV = _cached_junction_var_values[VolumeJunction1Phase::RHOEV_INDEX];
     145             : 
     146       69679 :   std::vector<ADReal> UJi(THMVACE3D::N_FLUX_INPUTS, 0.);
     147      139358 :   UJi[THMVACE3D::RHOA] = rhoV / _volume * Ai;
     148       69679 :   if (_apply_velocity_scaling)
     149             :   {
     150             :     const auto & rhoAi = Ui[THMVACE3D::RHOA];
     151             :     const auto & rhouAi = Ui[THMVACE3D::RHOUA];
     152             :     const auto & rhovAi = Ui[THMVACE3D::RHOVA];
     153             :     const auto & rhowAi = Ui[THMVACE3D::RHOWA];
     154             :     const auto & rhoEAi = Ui[THMVACE3D::RHOEA];
     155             : 
     156           0 :     const ADRealVectorValue uveci(rhouAi / rhoAi, rhovAi / rhoAi, rhowAi / rhoAi);
     157         476 :     const ADReal uni = uveci * nJi;
     158             : 
     159         476 :     const ADReal rhoJ = rhoV / _volume;
     160         952 :     const ADReal vJ = 1.0 / rhoJ;
     161           0 :     const ADRealVectorValue uvecJ(rhouV / rhoV, rhovV / rhoV, rhowV / rhoV);
     162         952 :     const ADReal eJ = rhoEV / rhoV - 0.5 * uvecJ * uvecJ;
     163             : 
     164         476 :     const ADReal unJ = uvecJ * nJi;
     165         476 :     const ADReal ut1J = uvecJ * t1;
     166         476 :     const ADReal ut2J = uvecJ * t2;
     167             : 
     168             :     const ADReal rhoi = rhoAi / Ai;
     169         476 :     const ADReal vi = 1.0 / rhoi;
     170         476 :     const ADReal ei = rhoEAi / rhoAi - 0.5 * uni * uni;
     171             : 
     172         476 :     const ADReal ci = _fp.c_from_v_e(vi, ei);
     173         476 :     const ADReal cJ = _fp.c_from_v_e(vJ, eJ);
     174         476 :     const ADReal cmax = max(ci, cJ);
     175             : 
     176         476 :     const ADReal uni_sign = (uni > 0) - (uni < 0);
     177        2380 :     const ADReal factor = 0.5 * (1.0 - uni_sign) * min(abs(uni - unJ) / cmax, 1.0);
     178             : 
     179         476 :     const ADReal unJ_mod = uni - factor * (uni - unJ);
     180         476 :     const ADRealVectorValue uvecJ_mod = unJ_mod * nJi + ut1J * t1 + ut2J * t2;
     181         952 :     const ADReal EJ_mod = eJ + 0.5 * uvecJ_mod * uvecJ_mod;
     182             : 
     183         476 :     UJi[THMVACE3D::RHOUA] = rhoJ * uvecJ_mod(0) * Ai;
     184         476 :     UJi[THMVACE3D::RHOVA] = rhoJ * uvecJ_mod(1) * Ai;
     185         476 :     UJi[THMVACE3D::RHOWA] = rhoJ * uvecJ_mod(2) * Ai;
     186         476 :     UJi[THMVACE3D::RHOEA] = rhoJ * EJ_mod * Ai;
     187             :   }
     188             :   else
     189             :   {
     190      138406 :     UJi[THMVACE3D::RHOUA] = rhouV / _volume * Ai;
     191      138406 :     UJi[THMVACE3D::RHOVA] = rhovV / _volume * Ai;
     192      138406 :     UJi[THMVACE3D::RHOWA] = rhowV / _volume * Ai;
     193      138406 :     UJi[THMVACE3D::RHOEA] = rhoEV / _volume * Ai;
     194             :   }
     195       69679 :   UJi[THMVACE3D::AREA] = Ai;
     196             : 
     197             :   std::vector<ADReal> FL, FR;
     198       69679 :   numerical_flux.calcFlux(UJi, Ui, nJi, t1, t2, FL, FR);
     199             : 
     200       69679 :   return FL;
     201       69679 : }
     202             : 
     203             : std::vector<ADReal>
     204       69461 : ADVolumeJunction1PhaseUserObject::computeResidual(const std::vector<ADReal> & flux_3d,
     205             :                                                   const std::vector<ADReal> & Ui,
     206             :                                                   const RealVectorValue & ni,
     207             :                                                   bool is_primary_connection) const
     208             : {
     209             :   using std::abs;
     210             : 
     211             :   const RealVectorValue nJi = -ni;
     212             : 
     213             :   RealVectorValue t1, t2;
     214       69461 :   THM::computeOrthogonalDirections(nJi, t1, t2);
     215             : 
     216             :   const auto & Ai = Ui[THMVACE3D::AREA];
     217             : 
     218             :   const auto & rhoV = _cached_junction_var_values[VolumeJunction1Phase::RHOV_INDEX];
     219             :   const auto & rhouV = _cached_junction_var_values[VolumeJunction1Phase::RHOUV_INDEX];
     220             :   const auto & rhovV = _cached_junction_var_values[VolumeJunction1Phase::RHOVV_INDEX];
     221             :   const auto & rhowV = _cached_junction_var_values[VolumeJunction1Phase::RHOWV_INDEX];
     222             :   const auto & rhoEV = _cached_junction_var_values[VolumeJunction1Phase::RHOEV_INDEX];
     223             : 
     224       69461 :   const ADReal rhoJ = rhoV / _volume;
     225      138922 :   const ADReal vJ = 1.0 / rhoJ;
     226       69461 :   const ADRealVectorValue uvecJ(rhouV / rhoV, rhovV / rhoV, rhowV / rhoV);
     227      138922 :   const ADReal eJ = rhoEV / rhoV - 0.5 * uvecJ * uvecJ;
     228       69461 :   const ADReal pJ = _fp.p_from_v_e(vJ, eJ);
     229             : 
     230       69461 :   std::vector<ADReal> residual(_n_scalar_eq, 0.0);
     231             : 
     232       69461 :   if (is_primary_connection && std::abs(_K) > 1e-10)
     233             :   {
     234             :     const auto & rhoAi = Ui[THMVACE3D::RHOA];
     235             :     const auto & rhouAi = Ui[THMVACE3D::RHOUA];
     236             :     const auto & rhovAi = Ui[THMVACE3D::RHOVA];
     237             :     const auto & rhowAi = Ui[THMVACE3D::RHOWA];
     238             :     const auto & rhoEAi = Ui[THMVACE3D::RHOEA];
     239             : 
     240           0 :     const ADRealVectorValue uveci(rhouAi / rhoAi, rhovAi / rhoAi, rhowAi / rhoAi);
     241        2321 :     const ADReal uni = uveci * nJi;
     242             : 
     243        2321 :     const ADReal v_in = THM::v_from_rhoA_A(rhoAi, Ai);
     244             :     const ADReal rhouA2 = rhouAi * rhouAi;
     245        2321 :     const ADReal e_in = rhoEAi / rhoAi - 0.5 * rhouA2 / (rhoAi * rhoAi);
     246        2321 :     const ADReal p_in = _fp.p_from_v_e(v_in, e_in);
     247        2321 :     const ADReal s0_in = _fp.s_from_v_e(v_in, e_in);
     248        2321 :     const ADReal T_in = _fp.T_from_v_e(v_in, e_in);
     249        2321 :     const ADReal h_in = _fp.h_from_p_T(p_in, T_in);
     250             :     const ADReal velin2 = uni * uni;
     251        2321 :     const ADReal h0_in = h_in + 0.5 * velin2;
     252        2321 :     const ADReal p0_in = _fp.p_from_h_s(h0_in, s0_in);
     253             :     ADReal S_loss;
     254        2321 :     if (_A_ref == 0)
     255        4138 :       S_loss = _K * (p0_in - p_in) * Ai;
     256             :     else
     257         504 :       S_loss = _K * (p0_in - p_in) * _A_ref;
     258        2321 :     if (uni > 0) // flow is out of the junction
     259             :     {
     260           0 :       residual[VolumeJunction1Phase::RHOUV_INDEX] -= ni(0) * S_loss;
     261           0 :       residual[VolumeJunction1Phase::RHOVV_INDEX] -= ni(1) * S_loss;
     262           0 :       residual[VolumeJunction1Phase::RHOWV_INDEX] -= ni(2) * S_loss;
     263             :     }
     264             :     else
     265             :     {
     266        2321 :       residual[VolumeJunction1Phase::RHOUV_INDEX] += ni(0) * S_loss;
     267        2321 :       residual[VolumeJunction1Phase::RHOVV_INDEX] += ni(1) * S_loss;
     268        2321 :       residual[VolumeJunction1Phase::RHOWV_INDEX] += ni(2) * S_loss;
     269             :     }
     270        4642 :     residual[VolumeJunction1Phase::RHOEV_INDEX] += S_loss * abs(uni);
     271             :   }
     272             : 
     273       69461 :   const ADRealVectorValue flux_mom_n = flux_3d[THMVACE3D::MOM_NORM] * nJi +
     274       69461 :                                        flux_3d[THMVACE3D::MOM_TAN1] * t1 +
     275       69461 :                                        flux_3d[THMVACE3D::MOM_TAN2] * t2;
     276             :   const RealVectorValue ex(1, 0, 0);
     277             :   const RealVectorValue ey(0, 1, 0);
     278             :   const RealVectorValue ez(0, 0, 1);
     279             : 
     280       69461 :   residual[VolumeJunction1Phase::RHOV_INDEX] -= -flux_3d[THMVACE3D::RHOA];
     281      138922 :   residual[VolumeJunction1Phase::RHOUV_INDEX] -= -flux_mom_n * ex - pJ * ni(0) * Ai;
     282      138922 :   residual[VolumeJunction1Phase::RHOVV_INDEX] -= -flux_mom_n * ey - pJ * ni(1) * Ai;
     283      138922 :   residual[VolumeJunction1Phase::RHOWV_INDEX] -= -flux_mom_n * ez - pJ * ni(2) * Ai;
     284       69461 :   residual[VolumeJunction1Phase::RHOEV_INDEX] -= -flux_3d[THMVACE3D::RHOEA];
     285             : 
     286       69461 :   return residual;
     287           0 : }
     288             : 
     289             : void
     290       35873 : ADVolumeJunction1PhaseUserObject::finalize()
     291             : {
     292       35873 :   ADVolumeJunctionBaseUserObject::finalize();
     293      215238 :   for (unsigned int i = 0; i < _n_scalar_eq; i++)
     294      179365 :     comm().sum(_residual[i]);
     295       35873 : }
     296             : 
     297             : std::vector<const MooseVariableBase *>
     298           0 : ADVolumeJunction1PhaseUserObject::getFlowChannelVariables() const
     299             : {
     300           0 :   std::vector<const MooseVariableBase *> vars(THMVACE1D::N_FLUX_OUTPUTS);
     301           0 :   vars[THMVACE1D::RHOA] = getVar("rhoA", 0);
     302           0 :   vars[THMVACE1D::RHOUA] = getVar("rhouA", 0);
     303           0 :   vars[THMVACE1D::RHOEA] = getVar("rhoEA", 0);
     304             : 
     305           0 :   return vars;
     306           0 : }
     307             : 
     308             : std::vector<const MooseVariableBase *>
     309       46056 : ADVolumeJunction1PhaseUserObject::getJunctionVariables() const
     310             : {
     311       46056 :   std::vector<const MooseVariableBase *> vars(VolumeJunction1Phase::N_EQ);
     312       46056 :   vars[VolumeJunction1Phase::RHOV_INDEX] = getJunctionVar("rhoV");
     313       46056 :   vars[VolumeJunction1Phase::RHOUV_INDEX] = getJunctionVar("rhouV");
     314       46056 :   vars[VolumeJunction1Phase::RHOVV_INDEX] = getJunctionVar("rhovV");
     315       46056 :   vars[VolumeJunction1Phase::RHOWV_INDEX] = getJunctionVar("rhowV");
     316       46056 :   vars[VolumeJunction1Phase::RHOEV_INDEX] = getJunctionVar("rhoEV");
     317             : 
     318       46056 :   return vars;
     319           0 : }