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 "VolumeJunction1Phase.h"
      11             : #include "FlowModelSinglePhase.h"
      12             : #include "THMMesh.h"
      13             : 
      14             : registerMooseObject("ThermalHydraulicsApp", VolumeJunction1Phase);
      15             : 
      16             : const unsigned int VolumeJunction1Phase::N_EQ = 5;
      17             : 
      18             : InputParameters
      19        1174 : VolumeJunction1Phase::validParams()
      20             : {
      21        1174 :   InputParameters params = FlowJunction1Phase::validParams();
      22             : 
      23        2348 :   params.addRequiredParam<Real>("volume", "Volume of the junction [m^3]");
      24        2348 :   params.addRequiredParam<Point>("position", "Spatial position of the center of the junction [m]");
      25             : 
      26        2348 :   params.addParam<FunctionName>("initial_p", "Initial pressure [Pa]");
      27        2348 :   params.addParam<FunctionName>("initial_T", "Initial temperature [K]");
      28        2348 :   params.addParam<FunctionName>("initial_vel_x", "Initial velocity in x-direction [m/s]");
      29        2348 :   params.addParam<FunctionName>("initial_vel_y", "Initial velocity in y-direction [m/s]");
      30        2348 :   params.addParam<FunctionName>("initial_vel_z", "Initial velocity in z-direction [m/s]");
      31             : 
      32        2348 :   params.addParam<Real>("scaling_factor_rhoV", 1.0, "Scaling factor for rho*V [-]");
      33        2348 :   params.addParam<Real>("scaling_factor_rhouV", 1.0, "Scaling factor for rho*u*V [-]");
      34        2348 :   params.addParam<Real>("scaling_factor_rhovV", 1.0, "Scaling factor for rho*v*V [-]");
      35        2348 :   params.addParam<Real>("scaling_factor_rhowV", 1.0, "Scaling factor for rho*w*V [-]");
      36        2348 :   params.addParam<Real>("scaling_factor_rhoEV", 1.0, "Scaling factor for rho*E*V [-]");
      37             : 
      38        2348 :   params.addParam<Real>("K", 0., "Form loss factor [-]");
      39        2348 :   params.addParam<Real>("A_ref", "Reference area [m^2]");
      40             : 
      41        3522 :   params.addParam<bool>("apply_velocity_scaling",
      42        2348 :                         false,
      43             :                         "Set to true to apply the scaling to the normal velocity. See "
      44             :                         "documentation for more information.");
      45             : 
      46        2348 :   params.declareControllable("K");
      47        1174 :   params.addClassDescription("Junction between 1-phase flow channels that has a non-zero volume");
      48             : 
      49        1174 :   return params;
      50           0 : }
      51             : 
      52         587 : VolumeJunction1Phase::VolumeJunction1Phase(const InputParameters & params)
      53             :   : FlowJunction1Phase(params),
      54             : 
      55         587 :     _volume(getParam<Real>("volume")),
      56        1174 :     _position(getParam<Point>("position")),
      57             : 
      58        1174 :     _scaling_factor_rhoV(getParam<Real>("scaling_factor_rhoV")),
      59        1174 :     _scaling_factor_rhouV(getParam<Real>("scaling_factor_rhouV")),
      60        1174 :     _scaling_factor_rhovV(getParam<Real>("scaling_factor_rhovV")),
      61        1174 :     _scaling_factor_rhowV(getParam<Real>("scaling_factor_rhowV")),
      62        1174 :     _scaling_factor_rhoEV(getParam<Real>("scaling_factor_rhoEV")),
      63             : 
      64         587 :     _rhoV_var_name(junctionVariableName("rhoV")),
      65         587 :     _rhouV_var_name(junctionVariableName("rhouV")),
      66         587 :     _rhovV_var_name(junctionVariableName("rhovV")),
      67         587 :     _rhowV_var_name(junctionVariableName("rhowV")),
      68         587 :     _rhoEV_var_name(junctionVariableName("rhoEV")),
      69         587 :     _pressure_var_name(junctionVariableName("p")),
      70         587 :     _temperature_var_name(junctionVariableName("T")),
      71         587 :     _velocity_var_name(junctionVariableName("vel")),
      72             : 
      73        1174 :     _K(getParam<Real>("K")),
      74        2149 :     _A_ref(isParamValid("A_ref") ? getParam<Real>("A_ref") : _zero)
      75             : {
      76             :   // Note: 'A_ref' can be required by child classes
      77        1200 :   if (!params.isParamRequired("A_ref") && params.isParamSetByUser("A_ref") &&
      78         639 :       !params.isParamSetByUser("K"))
      79           2 :     logWarning("Parameter 'A_ref' is specified, but 'K' is not specified, so the junction will "
      80             :                "behave as if there were no form loss.");
      81         587 : }
      82             : 
      83             : void
      84         587 : VolumeJunction1Phase::setupMesh()
      85             : {
      86         587 :   FlowJunction1Phase::setupMesh();
      87             : 
      88             :   // Add a NodeElem to the mesh
      89         587 :   auto * node = addNode(_position);
      90         587 :   auto * elem = addNodeElement(node->id());
      91         587 :   _junction_subdomain_id = mesh().getNextSubdomainId();
      92         587 :   elem->subdomain_id() = _junction_subdomain_id;
      93         587 :   setSubdomainInfo(_junction_subdomain_id, name());
      94             : 
      95             :   // Add coupling between the flow channel end elements and the NodeElem
      96             :   const auto & elem_ids = getConnectedElementIDs();
      97        1885 :   for (unsigned int i = 0; i < elem_ids.size(); i++)
      98        1298 :     getTHMProblem().augmentSparsity(elem_ids[i], elem->id());
      99         587 : }
     100             : 
     101             : void
     102         567 : VolumeJunction1Phase::check() const
     103             : {
     104         567 :   FlowJunction1Phase::check();
     105             : 
     106             :   bool ics_set =
     107        1129 :       getTHMProblem().hasInitialConditionsFromFile() ||
     108        3859 :       (isParamValid("initial_p") && isParamValid("initial_T") && isParamValid("initial_vel_x") &&
     109        2751 :        isParamValid("initial_vel_y") && isParamValid("initial_vel_z"));
     110             : 
     111         567 :   if (!ics_set && !_app.isRestarting())
     112             :   {
     113             :     // create a list of the missing IC parameters
     114             :     const std::vector<std::string> ic_params{
     115          24 :         "initial_p", "initial_T", "initial_vel_x", "initial_vel_y", "initial_vel_z"};
     116           4 :     std::ostringstream oss;
     117          24 :     for (const auto & ic_param : ic_params)
     118          20 :       if (!isParamValid(ic_param))
     119          20 :         oss << " " << ic_param;
     120             : 
     121           4 :     logError("The following initial condition parameters are missing:", oss.str());
     122           4 :   }
     123             : 
     124             :   // https://github.com/idaholab/moose/issues/28670
     125         567 :   if (getTHMProblem().hasInitialConditionsFromFile() && libMesh::n_threads() > 1 &&
     126           0 :       _app.n_processors() > 1)
     127           0 :     mooseDocumentedError("moose",
     128             :                          28670,
     129             :                          "Using initial conditions from a file for VolumeJunction1Phase is "
     130             :                          "currently not tested for parallel threading.");
     131         575 : }
     132             : 
     133             : void
     134         564 : VolumeJunction1Phase::addVariables()
     135             : {
     136         564 :   addJunctionVariable(true, _rhoV_var_name, _scaling_factor_rhoV);
     137         564 :   addJunctionVariable(true, _rhouV_var_name, _scaling_factor_rhouV);
     138         564 :   addJunctionVariable(true, _rhovV_var_name, _scaling_factor_rhovV);
     139         564 :   addJunctionVariable(true, _rhowV_var_name, _scaling_factor_rhowV);
     140         564 :   addJunctionVariable(true, _rhoEV_var_name, _scaling_factor_rhoEV);
     141             : 
     142         564 :   addJunctionVariable(false, _pressure_var_name);
     143         564 :   addJunctionVariable(false, _temperature_var_name);
     144         564 :   addJunctionVariable(false, _velocity_var_name);
     145             : 
     146        3863 :   if (isParamValid("initial_p") && isParamValid("initial_T") && isParamValid("initial_vel_x") &&
     147        3299 :       isParamValid("initial_vel_y") && isParamValid("initial_vel_z"))
     148             :   {
     149         547 :     addVolumeJunctionIC(_rhoV_var_name, "rhoV");
     150         547 :     addVolumeJunctionIC(_rhouV_var_name, "rhouV");
     151         547 :     addVolumeJunctionIC(_rhovV_var_name, "rhovV");
     152         547 :     addVolumeJunctionIC(_rhowV_var_name, "rhowV");
     153         547 :     addVolumeJunctionIC(_rhoEV_var_name, "rhoEV");
     154             : 
     155         547 :     addVolumeJunctionIC(_pressure_var_name, "p");
     156         547 :     addVolumeJunctionIC(_temperature_var_name, "T");
     157        1094 :     addVolumeJunctionIC(_velocity_var_name, "vel");
     158             :   }
     159         564 : }
     160             : 
     161             : void
     162         298 : VolumeJunction1Phase::buildVolumeJunctionUserObject()
     163             : {
     164         298 :   ExecFlagEnum execute_on(MooseUtils::getDefaultExecFlagEnum());
     165        1490 :   execute_on = {EXEC_INITIAL, EXEC_LINEAR, EXEC_NONLINEAR};
     166             : 
     167             :   {
     168         298 :     const std::string class_name = "ADVolumeJunction1PhaseUserObject";
     169         298 :     InputParameters params = _factory.getValidParams(class_name);
     170         298 :     params.set<bool>("use_scalar_variables") = false;
     171         298 :     params.set<subdomain_id_type>("junction_subdomain_id") = _junction_subdomain_id;
     172         298 :     params.set<std::vector<BoundaryName>>("boundary") = _boundary_names;
     173         596 :     params.set<std::vector<Real>>("normals") = _normals;
     174         298 :     params.set<std::vector<processor_id_type>>("processor_ids") = getConnectedProcessorIDs();
     175         298 :     params.set<std::vector<UserObjectName>>("numerical_flux_names") = _numerical_flux_names;
     176         298 :     params.set<Real>("volume") = _volume;
     177         894 :     params.set<std::vector<VariableName>>("A") = {FlowModel::AREA};
     178         894 :     params.set<std::vector<VariableName>>("rhoA") = {FlowModelSinglePhase::RHOA};
     179         894 :     params.set<std::vector<VariableName>>("rhouA") = {FlowModelSinglePhase::RHOUA};
     180         894 :     params.set<std::vector<VariableName>>("rhoEA") = {FlowModelSinglePhase::RHOEA};
     181         894 :     params.set<std::vector<VariableName>>("rhoV") = {_rhoV_var_name};
     182         894 :     params.set<std::vector<VariableName>>("rhouV") = {_rhouV_var_name};
     183         894 :     params.set<std::vector<VariableName>>("rhovV") = {_rhovV_var_name};
     184         894 :     params.set<std::vector<VariableName>>("rhowV") = {_rhowV_var_name};
     185         894 :     params.set<std::vector<VariableName>>("rhoEV") = {_rhoEV_var_name};
     186         298 :     params.set<Real>("K") = _K;
     187         298 :     params.set<Real>("A_ref") = _A_ref;
     188         298 :     params.set<UserObjectName>("fp") = _fp_name;
     189         596 :     params.set<bool>("apply_velocity_scaling") = getParam<bool>("apply_velocity_scaling");
     190         298 :     params.set<ExecFlagEnum>("execute_on") = execute_on;
     191         298 :     getTHMProblem().addUserObject(class_name, _junction_uo_name, params);
     192         298 :     connectObject(params, _junction_uo_name, "K");
     193         298 :   }
     194         596 : }
     195             : 
     196             : void
     197         564 : VolumeJunction1Phase::addMooseObjects()
     198             : {
     199         564 :   buildVolumeJunctionUserObject();
     200             : 
     201             :   // Add BC to each of the connected flow channels
     202        1805 :   for (std::size_t i = 0; i < _boundary_names.size(); i++)
     203             :   {
     204             :     const std::vector<NonlinearVariableName> var_names = {
     205        4964 :         FlowModelSinglePhase::RHOA, FlowModelSinglePhase::RHOUA, FlowModelSinglePhase::RHOEA};
     206        4964 :     for (std::size_t j = 0; j < var_names.size(); j++)
     207             :     {
     208        3723 :       const std::string class_name = "ADVolumeJunction1PhaseBC";
     209        3723 :       InputParameters params = _factory.getValidParams(class_name);
     210       11169 :       params.set<std::vector<BoundaryName>>("boundary") = {_boundary_names[i]};
     211        3723 :       params.set<Real>("normal") = _normals[i];
     212        3723 :       params.set<NonlinearVariableName>("variable") = var_names[j];
     213        7446 :       params.set<UserObjectName>("volume_junction_uo") = _junction_uo_name;
     214        3723 :       params.set<unsigned int>("connection_index") = i;
     215       11169 :       params.set<std::vector<VariableName>>("A_elem") = {FlowModel::AREA};
     216       11169 :       params.set<std::vector<VariableName>>("A_linear") = {FlowModel::AREA_LINEAR};
     217       11169 :       params.set<std::vector<VariableName>>("rhoA") = {FlowModelSinglePhase::RHOA};
     218       11169 :       params.set<std::vector<VariableName>>("rhouA") = {FlowModelSinglePhase::RHOUA};
     219       11169 :       params.set<std::vector<VariableName>>("rhoEA") = {FlowModelSinglePhase::RHOEA};
     220        3723 :       params.set<bool>("implicit") = getTHMProblem().getImplicitTimeIntegrationFlag();
     221        7446 :       getTHMProblem().addBoundaryCondition(
     222        3723 :           class_name, genName(name(), i, var_names[j] + ":" + class_name), params);
     223        3723 :     }
     224        1241 :   }
     225             : 
     226             :   // Add kernels for the junction
     227         564 :   std::vector<NonlinearVariableName> var_names(N_EQ);
     228         564 :   var_names[RHOV_INDEX] = _rhoV_var_name;
     229         564 :   var_names[RHOUV_INDEX] = _rhouV_var_name;
     230         564 :   var_names[RHOVV_INDEX] = _rhovV_var_name;
     231         564 :   var_names[RHOWV_INDEX] = _rhowV_var_name;
     232         564 :   var_names[RHOEV_INDEX] = _rhoEV_var_name;
     233        3384 :   for (std::size_t i = 0; i < N_EQ; i++)
     234             :   {
     235             :     {
     236        2820 :       const std::string class_name = "ADTimeDerivative";
     237        2820 :       InputParameters params = _factory.getValidParams(class_name);
     238        2820 :       params.set<NonlinearVariableName>("variable") = var_names[i];
     239        2820 :       const std::string obj_name = genName(name(), var_names[i], "td");
     240        2820 :       params.set<std::vector<SubdomainName>>("block") = getSubdomainNames();
     241        2820 :       getTHMProblem().addKernel(class_name, obj_name, params);
     242        2820 :     }
     243             :     {
     244        2820 :       const std::string class_name = "ADVolumeJunctionAdvectionKernel";
     245        2820 :       InputParameters params = _factory.getValidParams(class_name);
     246        2820 :       params.set<NonlinearVariableName>("variable") = var_names[i];
     247        5640 :       params.set<UserObjectName>("volume_junction_uo") = _junction_uo_name;
     248        2820 :       params.set<unsigned int>("equation_index") = i;
     249        2820 :       const std::string obj_name = genName(name(), var_names[i], "vja_sk");
     250        2820 :       params.set<std::vector<SubdomainName>>("block") = getSubdomainNames();
     251        2820 :       getTHMProblem().addKernel(class_name, obj_name, params);
     252        2820 :     }
     253             :   }
     254             : 
     255             :   const std::vector<std::pair<std::string, VariableName>> quantities = {
     256         564 :       {"pressure", _pressure_var_name},
     257         564 :       {"temperature", _temperature_var_name},
     258        2256 :       {"speed", _velocity_var_name}};
     259        2256 :   for (const auto & quantity_and_name : quantities)
     260             :   {
     261        1692 :     const std::string class_name = "VolumeJunction1PhaseAux";
     262        1692 :     InputParameters params = _factory.getValidParams(class_name);
     263        3384 :     params.set<AuxVariableName>("variable") = quantity_and_name.second;
     264        1692 :     params.set<MooseEnum>("quantity") = quantity_and_name.first;
     265        1692 :     params.set<Real>("volume") = _volume;
     266        5076 :     params.set<std::vector<VariableName>>("rhoV") = {_rhoV_var_name};
     267        5076 :     params.set<std::vector<VariableName>>("rhouV") = {_rhouV_var_name};
     268        5076 :     params.set<std::vector<VariableName>>("rhovV") = {_rhovV_var_name};
     269        5076 :     params.set<std::vector<VariableName>>("rhowV") = {_rhowV_var_name};
     270        5076 :     params.set<std::vector<VariableName>>("rhoEV") = {_rhoEV_var_name};
     271        1692 :     params.set<UserObjectName>("fp") = _fp_name;
     272        3384 :     const std::string obj_name = genName(name(), quantity_and_name.first + "_aux");
     273        1692 :     params.set<std::vector<SubdomainName>>("block") = getSubdomainNames();
     274        1692 :     getTHMProblem().addAuxKernel(class_name, obj_name, params);
     275        1692 :   }
     276             : 
     277             :   // An error message results if there is any block without a material, so
     278             :   // until this restriction is removed, we must add a dummy material that
     279             :   // computes no material properties.
     280             :   {
     281         564 :     const std::string class_name = "GenericConstantMaterial";
     282         564 :     InputParameters params = _factory.getValidParams(class_name);
     283         564 :     params.set<std::vector<SubdomainName>>("block") = getSubdomainNames();
     284         564 :     params.set<std::vector<std::string>>("prop_names") = {};
     285         564 :     params.set<std::vector<Real>>("prop_values") = {};
     286        1128 :     getTHMProblem().addMaterial(class_name, genName(name(), "dummy_mat"), params);
     287         564 :   }
     288        1128 : }
     289             : 
     290             : std::string
     291        5233 : VolumeJunction1Phase::junctionVariableName(const std::string & var_base) const
     292             : {
     293        5233 :   return var_base;
     294             : }
     295             : 
     296             : void
     297        5019 : VolumeJunction1Phase::addJunctionVariable(bool is_nonlinear,
     298             :                                           const VariableName & var,
     299             :                                           Real scaling_factor)
     300             : {
     301        5019 :   auto & problem = getTHMProblem();
     302             : 
     303             :   const libMesh::FEType fe_type(CONSTANT, MONOMIAL);
     304        5019 :   const auto & subdomains = getSubdomainNames();
     305             : 
     306        5019 :   if (is_nonlinear)
     307        2820 :     problem.addSimVariable(is_nonlinear, var, fe_type, subdomains, scaling_factor);
     308             :   else
     309        2199 :     problem.addSimVariable(is_nonlinear, var, fe_type, subdomains);
     310        5019 : }
     311             : 
     312             : void
     313         475 : VolumeJunction1Phase::addJunctionIC(const VariableName & var, Real value)
     314             : {
     315         475 :   getTHMProblem().addConstantIC(var, value, getSubdomainNames());
     316         475 : }
     317             : 
     318             : void
     319        4376 : VolumeJunction1Phase::addVolumeJunctionIC(const VariableName & var, const std::string & quantity)
     320             : {
     321        4376 :   const std::string class_name = "VolumeJunction1PhaseIC";
     322        4376 :   InputParameters params = _factory.getValidParams(class_name);
     323        4376 :   params.set<std::vector<SubdomainName>>("block") = getSubdomainNames();
     324        4376 :   params.set<VariableName>("variable") = var;
     325        4376 :   params.set<MooseEnum>("quantity") = quantity;
     326        4376 :   params.applySpecificParameters(parameters(),
     327             :                                  {"initial_p",
     328             :                                   "initial_T",
     329             :                                   "initial_vel_x",
     330             :                                   "initial_vel_y",
     331             :                                   "initial_vel_z",
     332             :                                   "volume",
     333             :                                   "position"});
     334        4376 :   params.set<UserObjectName>("fluid_properties") = _fp_name;
     335        8752 :   getTHMProblem().addSimInitialCondition(class_name, genName(name(), var, "ic"), params);
     336        8752 : }