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 "SCMQuadInterWrapperMeshGenerator.h"
      11             : #include "QuadInterWrapperMesh.h"
      12             : #include "libmesh/edge_edge2.h"
      13             : #include <numeric>
      14             : 
      15             : registerMooseObject("SubChannelApp", SCMQuadInterWrapperMeshGenerator);
      16             : registerMooseObjectRenamed("SubChannelApp",
      17             :                            QuadInterWrapperMeshGenerator,
      18             :                            "06/30/2025 24:00",
      19             :                            SCMQuadInterWrapperMeshGenerator);
      20             : 
      21             : InputParameters
      22          84 : SCMQuadInterWrapperMeshGenerator::validParams()
      23             : {
      24          84 :   InputParameters params = MeshGenerator::validParams();
      25          84 :   params.addClassDescription("Creates a mesh for the inter-wrapper around square subassemblies");
      26         168 :   params.addRequiredParam<Real>("assembly_pitch", "Assembly Pitch [m]");
      27         168 :   params.addRequiredParam<Real>("assembly_side_x",
      28             :                                 "Outer side lengths of assembly in x [m] - including duct");
      29         168 :   params.addRequiredParam<Real>("assembly_side_y",
      30             :                                 "Outer side lengths of assembly in y [m] - including duct");
      31         168 :   params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
      32         168 :   params.addRequiredParam<Real>("heated_length", "Heated length [m]");
      33         168 :   params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");
      34         168 :   params.addParam<Real>("Kij", 0.5, "Lateral form loss coefficient [-]");
      35         168 :   params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
      36         168 :   params.addRequiredParam<unsigned int>("nx", "Number of assemblies in the x direction [-]");
      37         168 :   params.addRequiredParam<unsigned int>("ny", "Number of assemblies in the y direction [-]");
      38         168 :   params.addRequiredParam<Real>("side_bypass",
      39             :                                 "Extra size of the bypass for the side assemblies [m]");
      40         168 :   params.addParam<unsigned int>("block_id", 0, "Domain Index");
      41          84 :   return params;
      42           0 : }
      43             : 
      44          42 : SCMQuadInterWrapperMeshGenerator::SCMQuadInterWrapperMeshGenerator(const InputParameters & params)
      45             :   : MeshGenerator(params),
      46          42 :     _assembly_pitch(getParam<Real>("assembly_pitch")),
      47          84 :     _assembly_side_x(getParam<Real>("assembly_side_x")),
      48          84 :     _assembly_side_y(getParam<Real>("assembly_side_y")),
      49          84 :     _unheated_length_entry(getParam<Real>("unheated_length_entry")),
      50          84 :     _heated_length(getParam<Real>("heated_length")),
      51          84 :     _unheated_length_exit(getParam<Real>("unheated_length_exit")),
      52          84 :     _kij(getParam<Real>("Kij")),
      53          84 :     _n_cells(getParam<unsigned int>("n_cells")),
      54          84 :     _nx(getParam<unsigned int>("nx")),
      55          84 :     _ny(getParam<unsigned int>("ny")),
      56          84 :     _side_bypass_length(getParam<Real>("side_bypass")),
      57          42 :     _n_channels((_nx + 1) * (_ny + 1)),
      58          42 :     _n_gaps(_nx * (_ny + 1) + _ny * (_nx + 1)),
      59          42 :     _n_assemblies(_nx * _ny),
      60         126 :     _block_id(getParam<unsigned int>("block_id"))
      61             : {
      62             :   // Converting number of assemblies into number of inter-wrapper flow channels
      63          42 :   _nx += 1;
      64          42 :   _ny += 1;
      65             : 
      66          42 :   if (_nx < 2 && _ny < 2)
      67           0 :     mooseError(name(), ": The number of assemblies cannot be less than one in any direction. ");
      68             : 
      69             :   // Defining the total length from 3 axial sections
      70          42 :   Real L = _unheated_length_entry + _heated_length + _unheated_length_exit;
      71             : 
      72             :   // Defining the dz based in the total length and the specified number of axial cells
      73          42 :   Real dz = L / _n_cells;
      74         864 :   for (unsigned int i = 0; i < _n_cells + 1; i++)
      75         822 :     _z_grid.push_back(dz * i);
      76             : 
      77             :   // Defining the array for axial resistances
      78          42 :   _k_grid.resize(_n_channels, std::vector<Real>(_n_cells + 1));
      79             : 
      80             :   // Defining the size of the maps
      81          42 :   _gap_to_chan_map.resize(_n_gaps);
      82          42 :   _gapnodes.resize(_n_gaps);
      83          42 :   _chan_to_gap_map.resize(_n_channels);
      84          42 :   _chan_to_pin_map.resize(_n_channels);
      85          42 :   _pin_to_chan_map.resize(_n_assemblies);
      86          42 :   _sign_id_crossflow_map.resize(_n_channels);
      87          42 :   _gij_map.resize(_n_gaps);
      88             : 
      89             :   // Defining the signs for positive and negative flows
      90          42 :   Real positive_flow = 1.0;
      91          42 :   Real negative_flow = -1.0;
      92             : 
      93             :   // Defining the subchannel types
      94          42 :   _subch_type.resize(_n_channels);
      95         366 :   for (unsigned int iy = 0; iy < _ny; iy++)
      96             :   {
      97        3636 :     for (unsigned int ix = 0; ix < _nx; ix++)
      98             :     {
      99        3312 :       unsigned int i_ch = _nx * iy + ix;
     100        3270 :       bool is_corner = (ix == 0 && iy == 0) || (ix == _nx - 1 && iy == 0) ||
     101        6540 :                        (ix == 0 && iy == _ny - 1) || (ix == _nx - 1 && iy == _ny - 1);
     102        3312 :       bool is_edge = (ix == 0 || iy == 0 || ix == _nx - 1 || iy == _ny - 1);
     103             : 
     104        3312 :       if (is_corner)
     105         168 :         _subch_type[i_ch] = EChannelType::CORNER;
     106        3144 :       else if (is_edge)
     107         960 :         _subch_type[i_ch] = EChannelType::EDGE;
     108             :       else
     109        2184 :         _subch_type[i_ch] = EChannelType::CENTER;
     110             :     }
     111             :   }
     112             : 
     113             :   // Index the east-west gaps.
     114          42 :   unsigned int i_gap = 0;
     115         366 :   for (unsigned int iy = 0; iy < _ny; iy++)
     116             :   {
     117        3312 :     for (unsigned int ix = 0; ix < _nx - 1; ix++)
     118             :     {
     119        2988 :       unsigned int i_ch = _nx * iy + ix;
     120        2988 :       unsigned int j_ch = _nx * iy + (ix + 1);
     121        2988 :       _gap_to_chan_map[i_gap] = {i_ch, j_ch};
     122        2988 :       _chan_to_gap_map[i_ch].push_back(i_gap);
     123        2988 :       _chan_to_gap_map[j_ch].push_back(i_gap);
     124        2988 :       _sign_id_crossflow_map[i_ch].push_back(positive_flow);
     125        2988 :       _sign_id_crossflow_map[j_ch].push_back(negative_flow);
     126             : 
     127             :       // make a gap size map
     128        2988 :       if (iy == 0 || iy == _ny - 1)
     129         564 :         _gij_map[i_gap] = (_assembly_pitch - _assembly_side_x) / 2 + _side_bypass_length;
     130             :       else
     131        2424 :         _gij_map[i_gap] = (_assembly_pitch - _assembly_side_x);
     132        2988 :       ++i_gap;
     133             :     }
     134             :   }
     135             : 
     136             :   // Index the north-south gaps.
     137         324 :   for (unsigned int iy = 0; iy < _ny - 1; iy++)
     138             :   {
     139        3270 :     for (unsigned int ix = 0; ix < _nx; ix++)
     140             :     {
     141        2988 :       unsigned int i_ch = _nx * iy + ix;
     142        2988 :       unsigned int j_ch = _nx * (iy + 1) + ix;
     143        2988 :       _gap_to_chan_map[i_gap] = {i_ch, j_ch};
     144        2988 :       _chan_to_gap_map[i_ch].push_back(i_gap);
     145        2988 :       _chan_to_gap_map[j_ch].push_back(i_gap);
     146        2988 :       _sign_id_crossflow_map[i_ch].push_back(positive_flow);
     147        2988 :       _sign_id_crossflow_map[j_ch].push_back(negative_flow);
     148             : 
     149             :       // make a gap size map
     150        2988 :       if (ix == 0 || ix == _nx - 1)
     151         564 :         _gij_map[i_gap] = (_assembly_pitch - _assembly_side_y) / 2 + _side_bypass_length;
     152             :       else
     153        2424 :         _gij_map[i_gap] = (_assembly_pitch - _assembly_side_y);
     154        2988 :       ++i_gap;
     155             :     }
     156             :   }
     157             : 
     158             :   // Make assembly to channel map
     159         324 :   for (unsigned int iy = 0; iy < _ny - 1; iy++)
     160             :   {
     161        2988 :     for (unsigned int ix = 0; ix < _nx - 1; ix++)
     162             :     {
     163        2706 :       unsigned int i_pin = (_nx - 1) * iy + ix;
     164        2706 :       unsigned int i_chan_1 = _nx * iy + ix;
     165        2706 :       unsigned int i_chan_2 = _nx * (iy + 1) + ix;
     166        2706 :       unsigned int i_chan_3 = _nx * (iy + 1) + (ix + 1);
     167        2706 :       unsigned int i_chan_4 = _nx * iy + (ix + 1);
     168        2706 :       _pin_to_chan_map[i_pin].push_back(i_chan_1);
     169        2706 :       _pin_to_chan_map[i_pin].push_back(i_chan_2);
     170        2706 :       _pin_to_chan_map[i_pin].push_back(i_chan_3);
     171        2706 :       _pin_to_chan_map[i_pin].push_back(i_chan_4);
     172             :     }
     173             :   }
     174             : 
     175             :   // Make channel to assembly map
     176         366 :   for (unsigned int iy = 0; iy < _ny; iy++) // row
     177             :   {
     178        3636 :     for (unsigned int ix = 0; ix < _nx; ix++) // column
     179             :     {
     180        3312 :       unsigned int i_ch = _nx * iy + ix;
     181             :       // Corners contact 1/4 of one pin
     182        3312 :       if (iy == 0 && ix == 0)
     183             :       {
     184          42 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     185             :       }
     186        3270 :       else if (iy == _ny - 1 && ix == 0)
     187             :       {
     188          42 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     189             :       }
     190        3228 :       else if (iy == 0 && ix == _nx - 1)
     191             :       {
     192          42 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     193             :       }
     194        3186 :       else if (iy == _ny - 1 && ix == _nx - 1)
     195             :       {
     196          42 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
     197             :       }
     198             :       // Sides contact 1/4 of two pins
     199        3144 :       else if (iy == 0)
     200             :       {
     201         240 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     202         240 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     203             :       }
     204        2904 :       else if (iy == _ny - 1)
     205             :       {
     206         240 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     207         240 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
     208             :       }
     209        2664 :       else if (ix == 0)
     210             :       {
     211         240 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     212         240 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     213             :       }
     214        2424 :       else if (ix == _nx - 1)
     215             :       {
     216         240 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     217         240 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
     218             :       }
     219             :       // interior contacts 1/4 of 4 pins
     220             :       else
     221             :       {
     222        2184 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix);
     223        2184 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * iy + ix - 1);
     224        2184 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix);
     225        2184 :         _chan_to_pin_map[i_ch].push_back((_nx - 1) * (iy - 1) + ix - 1);
     226             :       }
     227             :     }
     228             :   }
     229             : 
     230             :   // Reduce reserved memory in the channel-to-gap map.
     231        3354 :   for (auto & gap : _chan_to_gap_map)
     232             :     gap.shrink_to_fit();
     233             : 
     234             :   // Reduce reserved memory in the channel-to-pin map.
     235        3354 :   for (auto & pin : _chan_to_pin_map)
     236             :     pin.shrink_to_fit();
     237             : 
     238             :   // Reduce reserved memory in the pin-to-channel map.
     239        2748 :   for (auto & pin : _pin_to_chan_map)
     240             :     pin.shrink_to_fit();
     241             : 
     242          42 :   _console << "Inter-wrapper quad mesh initialized" << std::endl;
     243          42 : }
     244             : 
     245             : std::unique_ptr<MeshBase>
     246          42 : SCMQuadInterWrapperMeshGenerator::generate()
     247             : {
     248          42 :   auto mesh_base = buildMeshBaseObject();
     249             :   BoundaryInfo & boundary_info = mesh_base->get_boundary_info();
     250          42 :   mesh_base->set_spatial_dimension(3);
     251          42 :   mesh_base->reserve_elem(_n_cells * _ny * _nx);
     252          42 :   mesh_base->reserve_nodes((_n_cells + 1) * _ny * _nx);
     253          42 :   _nodes.resize(_nx * _ny);
     254             :   // Add the points in the shape of a rectilinear grid.  The grid is regular
     255             :   // on the xy-plane with a spacing of `pitch` between points.  The grid along
     256             :   // z is irregular to account for Pin spacers.  Store pointers in the _nodes
     257             :   // array so we can keep track of which points are in which channels.
     258          42 :   Real offset_x = (_nx - 1) * _assembly_pitch / 2.0;
     259          42 :   Real offset_y = (_ny - 1) * _assembly_pitch / 2.0;
     260             :   unsigned int node_id = 0;
     261         366 :   for (unsigned int iy = 0; iy < _ny; iy++)
     262             :   {
     263        3636 :     for (unsigned int ix = 0; ix < _nx; ix++)
     264             :     {
     265        3312 :       int i_ch = _nx * iy + ix;
     266        3312 :       _nodes[i_ch].reserve(_n_cells);
     267      131904 :       for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
     268             :       {
     269      128592 :         _nodes[i_ch].push_back(mesh_base->add_point(
     270      128592 :             Point(_assembly_pitch * ix - offset_x, _assembly_pitch * iy - offset_y, _z_grid[iz]),
     271      128592 :             node_id++));
     272             :       }
     273             :     }
     274             :   }
     275             : 
     276             :   // Add the elements which in this case are 2-node edges that link each
     277             :   // subchannel's nodes vertically.
     278             :   unsigned int elem_id = 0;
     279         366 :   for (unsigned int iy = 0; iy < _ny; iy++)
     280             :   {
     281        3636 :     for (unsigned int ix = 0; ix < _nx; ix++)
     282             :     {
     283      128592 :       for (unsigned int iz = 0; iz < _n_cells; iz++)
     284             :       {
     285      125280 :         Elem * elem = new Edge2;
     286      125280 :         elem->subdomain_id() = _block_id;
     287      125280 :         elem->set_id(elem_id++);
     288      125280 :         elem = mesh_base->add_elem(elem);
     289      125280 :         const int indx1 = ((_n_cells + 1) * _nx) * iy + (_n_cells + 1) * ix + iz;
     290      125280 :         const int indx2 = ((_n_cells + 1) * _nx) * iy + (_n_cells + 1) * ix + (iz + 1);
     291      125280 :         elem->set_node(0, mesh_base->node_ptr(indx1));
     292      125280 :         elem->set_node(1, mesh_base->node_ptr(indx2));
     293             : 
     294      125280 :         if (iz == 0)
     295        3312 :           boundary_info.add_side(elem, 0, 0);
     296      125280 :         if (iz == _n_cells - 1)
     297        3312 :           boundary_info.add_side(elem, 1, 1);
     298             :       }
     299             :     }
     300             :   }
     301             : 
     302          42 :   boundary_info.sideset_name(0) = "inlet";
     303          42 :   boundary_info.sideset_name(1) = "outlet";
     304          42 :   boundary_info.nodeset_name(0) = "inlet";
     305          42 :   boundary_info.nodeset_name(1) = "outlet";
     306          42 :   mesh_base->subdomain_name(_block_id) = name();
     307          42 :   mesh_base->prepare_for_use();
     308             : 
     309             :   // move the meta data into QuadInterWrapperMesh
     310          42 :   auto & sch_mesh = static_cast<QuadInterWrapperMesh &>(*_mesh);
     311          42 :   sch_mesh._unheated_length_entry = _unheated_length_entry;
     312             : 
     313          42 :   sch_mesh._assembly_pitch = _assembly_pitch;
     314          42 :   sch_mesh._assembly_side_x = _assembly_side_x;
     315          42 :   sch_mesh._assembly_side_y = _assembly_side_y;
     316          42 :   sch_mesh._heated_length = _heated_length;
     317          42 :   sch_mesh._unheated_length_exit = _unheated_length_exit;
     318          42 :   sch_mesh._z_grid = _z_grid;
     319          42 :   sch_mesh._k_grid = _k_grid;
     320          42 :   sch_mesh._kij = _kij;
     321          42 :   sch_mesh._n_cells = _n_cells;
     322          42 :   sch_mesh._nx = _nx;
     323          42 :   sch_mesh._ny = _ny;
     324          42 :   sch_mesh._side_bypass_length = _side_bypass_length;
     325          42 :   sch_mesh._n_channels = _n_channels;
     326          42 :   sch_mesh._n_gaps = _n_gaps;
     327          42 :   sch_mesh._n_assemblies = _n_assemblies;
     328          42 :   sch_mesh._nodes = _nodes;
     329          42 :   sch_mesh._gapnodes = _gapnodes;
     330          42 :   sch_mesh._gap_to_chan_map = _gap_to_chan_map;
     331          42 :   sch_mesh._chan_to_gap_map = _chan_to_gap_map;
     332          42 :   sch_mesh._chan_to_pin_map = _chan_to_pin_map;
     333          42 :   sch_mesh._pin_to_chan_map = _pin_to_chan_map;
     334          42 :   sch_mesh._sign_id_crossflow_map = _sign_id_crossflow_map;
     335          42 :   sch_mesh._gij_map = _gij_map;
     336          42 :   sch_mesh._subch_type = _subch_type;
     337             : 
     338          42 :   _console << "Inter-wrapper quad mesh generated" << std::endl;
     339             : 
     340          42 :   return mesh_base;
     341           0 : }