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      136247 : ADVolumeJunction1PhaseUserObject::computeFluxesAndResiduals(const unsigned int & c)
      98             : {
      99      136247 :   const Real din = _normal[c];
     100      136247 :   const RealVectorValue di = _dir[0];
     101             :   const RealVectorValue ni = di * din;
     102             :   const RealVectorValue nJi = -ni;
     103      136247 :   const Real nJi_dot_di = -din;
     104             : 
     105             :   RealVectorValue t1, t2;
     106      136247 :   THM::computeOrthogonalDirections(nJi, t1, t2);
     107             : 
     108      136247 :   std::vector<ADReal> Ui(THMVACE3D::N_FLUX_INPUTS, 0.);
     109      136247 :   Ui[THMVACE3D::RHOA] = _rhoA[0];
     110      272494 :   Ui[THMVACE3D::RHOUA] = _rhouA[0] * di(0);
     111      136247 :   Ui[THMVACE3D::RHOVA] = _rhouA[0] * di(1);
     112      136247 :   Ui[THMVACE3D::RHOWA] = _rhouA[0] * di(2);
     113      136247 :   Ui[THMVACE3D::RHOEA] = _rhoEA[0];
     114      136247 :   Ui[THMVACE3D::AREA] = _A[0];
     115             : 
     116             :   const auto & rhoV = _cached_junction_var_values[VolumeJunction1Phase::RHOV_INDEX];
     117             :   const auto & rhouV = _cached_junction_var_values[VolumeJunction1Phase::RHOUV_INDEX];
     118             :   const auto & rhovV = _cached_junction_var_values[VolumeJunction1Phase::RHOVV_INDEX];
     119             :   const auto & rhowV = _cached_junction_var_values[VolumeJunction1Phase::RHOWV_INDEX];
     120             :   const auto & rhoEV = _cached_junction_var_values[VolumeJunction1Phase::RHOEV_INDEX];
     121             : 
     122      136247 :   const ADReal rhoJ = rhoV / _volume;
     123      272494 :   const ADReal vJ = 1.0 / rhoJ;
     124           0 :   const ADRealVectorValue uvecJ(rhouV / rhoV, rhovV / rhoV, rhowV / rhoV);
     125      272494 :   const ADReal eJ = rhoEV / rhoV - 0.5 * uvecJ * uvecJ;
     126      136247 :   const ADReal pJ = _fp.p_from_v_e(vJ, eJ);
     127             : 
     128      136247 :   std::vector<ADReal> UJi(THMVACE3D::N_FLUX_INPUTS, 0.);
     129      272494 :   UJi[THMVACE3D::RHOA] = rhoV / _volume * _A[0];
     130      136247 :   if (_apply_velocity_scaling)
     131             :   {
     132           0 :     const ADReal unJ = uvecJ * nJi;
     133           0 :     const ADReal ut1J = uvecJ * t1;
     134           0 :     const ADReal ut2J = uvecJ * t2;
     135           0 :     const ADReal uni = _rhouA[0] / _rhoA[0] * nJi_dot_di;
     136             : 
     137           0 :     const ADReal rhoi = _rhoA[0] / _A[0];
     138           0 :     const ADReal vi = 1.0 / rhoi;
     139           0 :     const ADReal ei = _rhoEA[0] / _rhoA[0] - 0.5 * uni * uni;
     140           0 :     const ADReal ci = _fp.c_from_v_e(vi, ei);
     141           0 :     const ADReal cJ = _fp.c_from_v_e(vJ, eJ);
     142           0 :     const ADReal cmax = std::max(ci, cJ);
     143             : 
     144           0 :     const ADReal uni_sign = (uni > 0) - (uni < 0);
     145           0 :     const ADReal factor = 0.5 * (1.0 - uni_sign) * std::min(std::abs(uni - unJ) / cmax, 1.0);
     146             : 
     147           0 :     const ADReal unJ_mod = uni - factor * (uni - unJ);
     148           0 :     const ADRealVectorValue uvecJ_mod = unJ_mod * nJi + ut1J * t1 + ut2J * t2;
     149           0 :     const ADReal EJ_mod = eJ + 0.5 * uvecJ_mod * uvecJ_mod;
     150             : 
     151           0 :     UJi[THMVACE3D::RHOUA] = rhoJ * uvecJ_mod(0) * _A[0];
     152           0 :     UJi[THMVACE3D::RHOVA] = rhoJ * uvecJ_mod(1) * _A[0];
     153           0 :     UJi[THMVACE3D::RHOWA] = rhoJ * uvecJ_mod(2) * _A[0];
     154           0 :     UJi[THMVACE3D::RHOEA] = rhoJ * EJ_mod * _A[0];
     155             :   }
     156             :   else
     157             :   {
     158      272494 :     UJi[THMVACE3D::RHOUA] = rhouV / _volume * _A[0];
     159      272494 :     UJi[THMVACE3D::RHOVA] = rhovV / _volume * _A[0];
     160      272494 :     UJi[THMVACE3D::RHOWA] = rhowV / _volume * _A[0];
     161      272494 :     UJi[THMVACE3D::RHOEA] = rhoEV / _volume * _A[0];
     162             :   }
     163      136247 :   UJi[THMVACE3D::AREA] = _A[0];
     164             : 
     165             :   const auto flux_3d =
     166      136247 :       _numerical_flux_uo[c]->getFlux3D(_current_side, _current_elem->id(), UJi, Ui, nJi, t1, t2);
     167             : 
     168      136247 :   _flux[c].resize(THMVACE1D::N_FLUX_OUTPUTS);
     169      136247 :   _flux[c][THMVACE1D::MASS] = flux_3d[THMVACE3D::MASS] * nJi_dot_di;
     170      136247 :   _flux[c][THMVACE1D::MOMENTUM] = flux_3d[THMVACE3D::MOM_NORM];
     171      136247 :   _flux[c][THMVACE1D::ENERGY] = flux_3d[THMVACE3D::ENERGY] * nJi_dot_di;
     172             : 
     173      136247 :   if (c == 0 && std::abs(_K) > 1e-10)
     174             :   {
     175        3798 :     const ADReal vel_in = _rhouA[0] / _rhoA[0];
     176        3798 :     const ADReal v_in = THM::v_from_rhoA_A(_rhoA[0], _A[0]);
     177        3798 :     const ADReal rhouA2 = _rhouA[0] * _rhouA[0];
     178       11394 :     const ADReal e_in = _rhoEA[0] / _rhoA[0] - 0.5 * rhouA2 / (_rhoA[0] * _rhoA[0]);
     179        3798 :     const ADReal p_in = _fp.p_from_v_e(v_in, e_in);
     180        3798 :     const ADReal s0_in = _fp.s_from_v_e(v_in, e_in);
     181        3798 :     const ADReal T_in = _fp.T_from_v_e(v_in, e_in);
     182        3798 :     const ADReal h_in = _fp.h_from_p_T(p_in, T_in);
     183             :     const ADReal velin2 = vel_in * vel_in;
     184        3798 :     const ADReal h0_in = h_in + 0.5 * velin2;
     185        3798 :     const ADReal p0_in = _fp.p_from_h_s(h0_in, s0_in);
     186             :     ADReal S_loss;
     187        3798 :     if (_A_ref == 0)
     188       10449 :       S_loss = _K * (p0_in - p_in) * _A[0];
     189             :     else
     190         630 :       S_loss = _K * (p0_in - p_in) * _A_ref;
     191        3798 :     if (THM::isInlet(vel_in, _normal[c]))
     192             :     {
     193           0 :       _residual[VolumeJunction1Phase::RHOUV_INDEX] -= ni(0) * S_loss;
     194           0 :       _residual[VolumeJunction1Phase::RHOVV_INDEX] -= ni(1) * S_loss;
     195           0 :       _residual[VolumeJunction1Phase::RHOWV_INDEX] -= ni(2) * S_loss;
     196             :     }
     197             :     else
     198             :     {
     199        3798 :       _residual[VolumeJunction1Phase::RHOUV_INDEX] += ni(0) * S_loss;
     200        3798 :       _residual[VolumeJunction1Phase::RHOVV_INDEX] += ni(1) * S_loss;
     201        3798 :       _residual[VolumeJunction1Phase::RHOWV_INDEX] += ni(2) * S_loss;
     202             :     }
     203        7596 :     _residual[VolumeJunction1Phase::RHOEV_INDEX] += S_loss * std::abs(vel_in);
     204             :   }
     205             : 
     206      136247 :   const ADRealVectorValue flux_mom_n = flux_3d[THMVACE3D::MOM_NORM] * nJi +
     207      136247 :                                        flux_3d[THMVACE3D::MOM_TAN1] * t1 +
     208      136247 :                                        flux_3d[THMVACE3D::MOM_TAN2] * t2;
     209             :   const RealVectorValue ex(1, 0, 0);
     210             :   const RealVectorValue ey(0, 1, 0);
     211             :   const RealVectorValue ez(0, 0, 1);
     212             : 
     213      136247 :   _residual[VolumeJunction1Phase::RHOV_INDEX] -= -flux_3d[THMVACE3D::RHOA];
     214      408741 :   _residual[VolumeJunction1Phase::RHOUV_INDEX] -= -flux_mom_n * ex - pJ * ni(0) * _A[0];
     215      408741 :   _residual[VolumeJunction1Phase::RHOVV_INDEX] -= -flux_mom_n * ey - pJ * ni(1) * _A[0];
     216      408741 :   _residual[VolumeJunction1Phase::RHOWV_INDEX] -= -flux_mom_n * ez - pJ * ni(2) * _A[0];
     217      136247 :   _residual[VolumeJunction1Phase::RHOEV_INDEX] -= -flux_3d[THMVACE3D::RHOEA];
     218      136247 : }
     219             : 
     220             : void
     221       73760 : ADVolumeJunction1PhaseUserObject::finalize()
     222             : {
     223       73760 :   ADVolumeJunctionBaseUserObject::finalize();
     224      442560 :   for (unsigned int i = 0; i < _n_scalar_eq; i++)
     225      368800 :     comm().sum(_residual[i]);
     226       73760 : }