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        2372 : ADVolumeJunction1PhaseUserObject::validParams()
      26             : {
      27        2372 :   InputParameters params = ADVolumeJunctionBaseUserObject::validParams();
      28             : 
      29        4744 :   params.addRequiredCoupledVar("A", "Cross-sectional area of connected flow channels");
      30        4744 :   params.addRequiredCoupledVar("rhoA", "rho*A of the connected flow channels");
      31        4744 :   params.addRequiredCoupledVar("rhouA", "rhou*A of the connected flow channels");
      32        4744 :   params.addRequiredCoupledVar("rhoEA", "rhoE*A of the connected flow channels");
      33             : 
      34        4744 :   params.addRequiredCoupledVar("rhoV", "rho*V of the junction");
      35        4744 :   params.addRequiredCoupledVar("rhouV", "rho*u*V of the junction");
      36        4744 :   params.addRequiredCoupledVar("rhovV", "rho*v*V of the junction");
      37        4744 :   params.addRequiredCoupledVar("rhowV", "rho*w*V of the junction");
      38        4744 :   params.addRequiredCoupledVar("rhoEV", "rho*E*V of the junction");
      39             : 
      40        4744 :   params.addRequiredParam<UserObjectName>("fp", "Fluid properties user object name");
      41             : 
      42        4744 :   params.addRequiredParam<Real>("K", "Form loss coefficient [-]");
      43        4744 :   params.addRequiredParam<Real>("A_ref", "Reference area [m^2]");
      44             : 
      45        4744 :   params.addParam<bool>("apply_velocity_scaling",
      46        4744 :                         false,
      47             :                         "Set to true to apply the scaling to the normal velocity. See "
      48             :                         "documentation for more information.");
      49             : 
      50        4744 :   params.addClassDescription(
      51             :       "Computes and caches flux and residual vectors for a 1-phase volume junction");
      52             : 
      53        4744 :   params.declareControllable("K");
      54        2372 :   return params;
      55           0 : }
      56             : 
      57        1298 : ADVolumeJunction1PhaseUserObject::ADVolumeJunction1PhaseUserObject(const InputParameters & params)
      58             :   : ADVolumeJunctionBaseUserObject(params),
      59             : 
      60        1298 :     _A(adCoupledValue("A")),
      61        1298 :     _rhoA(adCoupledValue("rhoA")),
      62        1298 :     _rhouA(adCoupledValue("rhouA")),
      63        1298 :     _rhoEA(adCoupledValue("rhoEA")),
      64             : 
      65        2596 :     _K(getParam<Real>("K")),
      66        2596 :     _A_ref(getParam<Real>("A_ref")),
      67             : 
      68        2596 :     _apply_velocity_scaling(getParam<bool>("apply_velocity_scaling")),
      69             : 
      70        2596 :     _fp(getUserObject<SinglePhaseFluidProperties>("fp"))
      71             : {
      72        1298 :   _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        1298 :   _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        1298 :   _junction_var_values.resize(VolumeJunction1Phase::N_EQ);
      85        1298 :   _junction_var_values[VolumeJunction1Phase::RHOV_INDEX] = &coupledJunctionValue("rhoV");
      86        1298 :   _junction_var_values[VolumeJunction1Phase::RHOUV_INDEX] = &coupledJunctionValue("rhouV");
      87        1298 :   _junction_var_values[VolumeJunction1Phase::RHOVV_INDEX] = &coupledJunctionValue("rhovV");
      88        1298 :   _junction_var_values[VolumeJunction1Phase::RHOWV_INDEX] = &coupledJunctionValue("rhowV");
      89        1298 :   _junction_var_values[VolumeJunction1Phase::RHOEV_INDEX] = &coupledJunctionValue("rhoEV");
      90             : 
      91        1298 :   _numerical_flux_uo.resize(_n_connections);
      92        4174 :   for (std::size_t i = 0; i < _n_connections; i++)
      93        2876 :     _numerical_flux_uo[i] = &getUserObjectByName<ADNumericalFlux3EqnBase>(_numerical_flux_names[i]);
      94        1298 : }
      95             : 
      96             : void
      97      134457 : ADVolumeJunction1PhaseUserObject::computeFluxesAndResiduals(const unsigned int & c)
      98             : {
      99             :   using std::abs, std::max, std::min;
     100             : 
     101      134457 :   const Real din = _normal[c];
     102      134457 :   const RealVectorValue di = _dir[0];
     103             :   const RealVectorValue ni = di * din;
     104             :   const RealVectorValue nJi = -ni;
     105      134457 :   const Real nJi_dot_di = -din;
     106             : 
     107             :   RealVectorValue t1, t2;
     108      134457 :   THM::computeOrthogonalDirections(nJi, t1, t2);
     109             : 
     110      134457 :   std::vector<ADReal> Ui(THMVACE3D::N_FLUX_INPUTS, 0.);
     111      134457 :   Ui[THMVACE3D::RHOA] = _rhoA[0];
     112      268914 :   Ui[THMVACE3D::RHOUA] = _rhouA[0] * di(0);
     113      134457 :   Ui[THMVACE3D::RHOVA] = _rhouA[0] * di(1);
     114      134457 :   Ui[THMVACE3D::RHOWA] = _rhouA[0] * di(2);
     115      134457 :   Ui[THMVACE3D::RHOEA] = _rhoEA[0];
     116      134457 :   Ui[THMVACE3D::AREA] = _A[0];
     117             : 
     118             :   const auto & rhoV = _cached_junction_var_values[VolumeJunction1Phase::RHOV_INDEX];
     119             :   const auto & rhouV = _cached_junction_var_values[VolumeJunction1Phase::RHOUV_INDEX];
     120             :   const auto & rhovV = _cached_junction_var_values[VolumeJunction1Phase::RHOVV_INDEX];
     121             :   const auto & rhowV = _cached_junction_var_values[VolumeJunction1Phase::RHOWV_INDEX];
     122             :   const auto & rhoEV = _cached_junction_var_values[VolumeJunction1Phase::RHOEV_INDEX];
     123             : 
     124      134457 :   const ADReal rhoJ = rhoV / _volume;
     125      268914 :   const ADReal vJ = 1.0 / rhoJ;
     126           0 :   const ADRealVectorValue uvecJ(rhouV / rhoV, rhovV / rhoV, rhowV / rhoV);
     127      268914 :   const ADReal eJ = rhoEV / rhoV - 0.5 * uvecJ * uvecJ;
     128      134457 :   const ADReal pJ = _fp.p_from_v_e(vJ, eJ);
     129             : 
     130      134457 :   std::vector<ADReal> UJi(THMVACE3D::N_FLUX_INPUTS, 0.);
     131      268914 :   UJi[THMVACE3D::RHOA] = rhoV / _volume * _A[0];
     132      134457 :   if (_apply_velocity_scaling)
     133             :   {
     134           0 :     const ADReal unJ = uvecJ * nJi;
     135           0 :     const ADReal ut1J = uvecJ * t1;
     136           0 :     const ADReal ut2J = uvecJ * t2;
     137           0 :     const ADReal uni = _rhouA[0] / _rhoA[0] * nJi_dot_di;
     138             : 
     139           0 :     const ADReal rhoi = _rhoA[0] / _A[0];
     140           0 :     const ADReal vi = 1.0 / rhoi;
     141           0 :     const ADReal ei = _rhoEA[0] / _rhoA[0] - 0.5 * uni * uni;
     142           0 :     const ADReal ci = _fp.c_from_v_e(vi, ei);
     143           0 :     const ADReal cJ = _fp.c_from_v_e(vJ, eJ);
     144           0 :     const ADReal cmax = max(ci, cJ);
     145             : 
     146           0 :     const ADReal uni_sign = (uni > 0) - (uni < 0);
     147           0 :     const ADReal factor = 0.5 * (1.0 - uni_sign) * min(abs(uni - unJ) / cmax, 1.0);
     148             : 
     149           0 :     const ADReal unJ_mod = uni - factor * (uni - unJ);
     150           0 :     const ADRealVectorValue uvecJ_mod = unJ_mod * nJi + ut1J * t1 + ut2J * t2;
     151           0 :     const ADReal EJ_mod = eJ + 0.5 * uvecJ_mod * uvecJ_mod;
     152             : 
     153           0 :     UJi[THMVACE3D::RHOUA] = rhoJ * uvecJ_mod(0) * _A[0];
     154           0 :     UJi[THMVACE3D::RHOVA] = rhoJ * uvecJ_mod(1) * _A[0];
     155           0 :     UJi[THMVACE3D::RHOWA] = rhoJ * uvecJ_mod(2) * _A[0];
     156           0 :     UJi[THMVACE3D::RHOEA] = rhoJ * EJ_mod * _A[0];
     157             :   }
     158             :   else
     159             :   {
     160      268914 :     UJi[THMVACE3D::RHOUA] = rhouV / _volume * _A[0];
     161      268914 :     UJi[THMVACE3D::RHOVA] = rhovV / _volume * _A[0];
     162      268914 :     UJi[THMVACE3D::RHOWA] = rhowV / _volume * _A[0];
     163      268914 :     UJi[THMVACE3D::RHOEA] = rhoEV / _volume * _A[0];
     164             :   }
     165      134457 :   UJi[THMVACE3D::AREA] = _A[0];
     166             : 
     167             :   const auto flux_3d =
     168      134457 :       _numerical_flux_uo[c]->getFlux3D(_current_side, _current_elem->id(), UJi, Ui, nJi, t1, t2);
     169             : 
     170      134457 :   _flux[c].resize(THMVACE1D::N_FLUX_OUTPUTS);
     171      134457 :   _flux[c][THMVACE1D::MASS] = flux_3d[THMVACE3D::MASS] * nJi_dot_di;
     172      134457 :   _flux[c][THMVACE1D::MOMENTUM] = flux_3d[THMVACE3D::MOM_NORM];
     173      134457 :   _flux[c][THMVACE1D::ENERGY] = flux_3d[THMVACE3D::ENERGY] * nJi_dot_di;
     174             : 
     175      134457 :   if (c == 0 && abs(_K) > 1e-10)
     176             :   {
     177        3768 :     const ADReal vel_in = _rhouA[0] / _rhoA[0];
     178        3768 :     const ADReal v_in = THM::v_from_rhoA_A(_rhoA[0], _A[0]);
     179        3768 :     const ADReal rhouA2 = _rhouA[0] * _rhouA[0];
     180       11304 :     const ADReal e_in = _rhoEA[0] / _rhoA[0] - 0.5 * rhouA2 / (_rhoA[0] * _rhoA[0]);
     181        3768 :     const ADReal p_in = _fp.p_from_v_e(v_in, e_in);
     182        3768 :     const ADReal s0_in = _fp.s_from_v_e(v_in, e_in);
     183        3768 :     const ADReal T_in = _fp.T_from_v_e(v_in, e_in);
     184        3768 :     const ADReal h_in = _fp.h_from_p_T(p_in, T_in);
     185             :     const ADReal velin2 = vel_in * vel_in;
     186        3768 :     const ADReal h0_in = h_in + 0.5 * velin2;
     187        3768 :     const ADReal p0_in = _fp.p_from_h_s(h0_in, s0_in);
     188             :     ADReal S_loss;
     189        3768 :     if (_A_ref == 0)
     190       10359 :       S_loss = _K * (p0_in - p_in) * _A[0];
     191             :     else
     192         630 :       S_loss = _K * (p0_in - p_in) * _A_ref;
     193        3768 :     if (THM::isInlet(vel_in, _normal[c]))
     194             :     {
     195           0 :       _residual[VolumeJunction1Phase::RHOUV_INDEX] -= ni(0) * S_loss;
     196           0 :       _residual[VolumeJunction1Phase::RHOVV_INDEX] -= ni(1) * S_loss;
     197           0 :       _residual[VolumeJunction1Phase::RHOWV_INDEX] -= ni(2) * S_loss;
     198             :     }
     199             :     else
     200             :     {
     201        3768 :       _residual[VolumeJunction1Phase::RHOUV_INDEX] += ni(0) * S_loss;
     202        3768 :       _residual[VolumeJunction1Phase::RHOVV_INDEX] += ni(1) * S_loss;
     203        3768 :       _residual[VolumeJunction1Phase::RHOWV_INDEX] += ni(2) * S_loss;
     204             :     }
     205        7536 :     _residual[VolumeJunction1Phase::RHOEV_INDEX] += S_loss * abs(vel_in);
     206             :   }
     207             : 
     208      134457 :   const ADRealVectorValue flux_mom_n = flux_3d[THMVACE3D::MOM_NORM] * nJi +
     209      134457 :                                        flux_3d[THMVACE3D::MOM_TAN1] * t1 +
     210      134457 :                                        flux_3d[THMVACE3D::MOM_TAN2] * t2;
     211             :   const RealVectorValue ex(1, 0, 0);
     212             :   const RealVectorValue ey(0, 1, 0);
     213             :   const RealVectorValue ez(0, 0, 1);
     214             : 
     215      134457 :   _residual[VolumeJunction1Phase::RHOV_INDEX] -= -flux_3d[THMVACE3D::RHOA];
     216      403371 :   _residual[VolumeJunction1Phase::RHOUV_INDEX] -= -flux_mom_n * ex - pJ * ni(0) * _A[0];
     217      403371 :   _residual[VolumeJunction1Phase::RHOVV_INDEX] -= -flux_mom_n * ey - pJ * ni(1) * _A[0];
     218      403371 :   _residual[VolumeJunction1Phase::RHOWV_INDEX] -= -flux_mom_n * ez - pJ * ni(2) * _A[0];
     219      134457 :   _residual[VolumeJunction1Phase::RHOEV_INDEX] -= -flux_3d[THMVACE3D::RHOEA];
     220      134457 : }
     221             : 
     222             : void
     223       73114 : ADVolumeJunction1PhaseUserObject::finalize()
     224             : {
     225       73114 :   ADVolumeJunctionBaseUserObject::finalize();
     226      438684 :   for (unsigned int i = 0; i < _n_scalar_eq; i++)
     227      365570 :     comm().sum(_residual[i]);
     228       73114 : }