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 "SCMTriDuctMeshGenerator.h"
      11             : #include "TriSubChannelMesh.h"
      12             : #include <cmath>
      13             : #include "libmesh/face_quad4.h"
      14             : #include "libmesh/unstructured_mesh.h"
      15             : 
      16             : registerMooseObject("SubChannelApp", SCMTriDuctMeshGenerator);
      17             : 
      18             : InputParameters
      19         154 : SCMTriDuctMeshGenerator::validParams()
      20             : {
      21         154 :   InputParameters params = MeshGenerator::validParams();
      22         154 :   params.addClassDescription(
      23             :       "Creates a mesh of 2D duct cells around a triangular lattice subassembly");
      24         308 :   params.addRequiredParam<MeshGeneratorName>("input", "The corresponding subchannel mesh");
      25         308 :   params.addParam<unsigned int>("block_id", 2, "Domain Index");
      26         308 :   params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
      27         308 :   params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
      28         308 :   params.addRequiredParam<Real>("heated_length", "Heated length [m]");
      29         308 :   params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");
      30         308 :   params.addRequiredParam<Real>("pitch", "Pitch is the distance between adjacent pins [m]");
      31         308 :   params.addRequiredParam<unsigned int>("nrings", "Number of fuel Pin rings in the assembly [-]");
      32         308 :   params.addRequiredParam<Real>("flat_to_flat",
      33             :                                 "Flat to flat distance for the hexagonal assembly [m]");
      34         154 :   return params;
      35           0 : }
      36             : 
      37          77 : SCMTriDuctMeshGenerator::SCMTriDuctMeshGenerator(const InputParameters & params)
      38             :   : MeshGenerator(params),
      39          77 :     _input(getMesh("input")),
      40         154 :     _n_cells(getParam<unsigned int>("n_cells")),
      41         154 :     _unheated_length_entry(getParam<Real>("unheated_length_entry")),
      42         154 :     _heated_length(getParam<Real>("heated_length")),
      43         154 :     _unheated_length_exit(getParam<Real>("unheated_length_exit")),
      44         154 :     _block_id(getParam<unsigned int>("block_id")),
      45         154 :     _pitch(getParam<Real>("pitch")),
      46         154 :     _n_rings(getParam<unsigned int>("nrings")),
      47         231 :     _flat_to_flat(getParam<Real>("flat_to_flat"))
      48             : {
      49          77 :   SubChannelMesh::generateZGrid(
      50          77 :       _unheated_length_entry, _heated_length, _unheated_length_exit, _n_cells, _z_grid);
      51          77 : }
      52             : 
      53             : std::unique_ptr<MeshBase>
      54          77 : SCMTriDuctMeshGenerator::generate()
      55             : {
      56          77 :   std::unique_ptr<MeshBase> mesh_base = std::move(_input);
      57          77 :   mesh_base->set_mesh_dimension(3);
      58             : 
      59             :   std::vector<Point> corners;
      60          77 :   ductCorners(corners, _flat_to_flat, Point(0, 0, 0));
      61             :   std::vector<Point> cross_sec;
      62          77 :   ductCrossSec(cross_sec, corners, _n_rings, _pitch);
      63             :   std::vector<Point> points;
      64          77 :   ductPoints(points, cross_sec, _z_grid);
      65             :   std::vector<std::vector<size_t>> elem_point_indices;
      66          77 :   ductElems(elem_point_indices, _z_grid.size(), cross_sec.size());
      67             :   std::vector<Node *> duct_nodes;
      68          77 :   buildDuct(mesh_base, duct_nodes, points, elem_point_indices, _block_id);
      69          77 :   mesh_base->subdomain_name(_block_id) = name();
      70             : 
      71          77 :   mesh_base->prepare_for_use();
      72             : 
      73          77 :   auto & sch_mesh = static_cast<TriSubChannelMesh &>(*_mesh);
      74          77 :   sch_mesh.setChannelToDuctMaps(duct_nodes);
      75          77 :   sch_mesh._duct_mesh_exist = true;
      76             : 
      77          77 :   return mesh_base;
      78          77 : }
      79             : 
      80             : void
      81          77 : SCMTriDuctMeshGenerator::ductCorners(std::vector<Point> & corners,
      82             :                                      Real flat_to_flat,
      83             :                                      const Point & center) const
      84             : {
      85          77 :   corners.resize(TriSubChannelMesh::N_CORNERS);
      86          77 :   Real r_corner = flat_to_flat / 2. / std::cos(libMesh::pi / 6.);
      87         539 :   for (size_t i = 0; i < TriSubChannelMesh::N_CORNERS; i++)
      88             :   {
      89         462 :     Real theta = i * libMesh::pi / 3.;
      90         462 :     Point corner = {r_corner * std::cos(theta), r_corner * std::sin(theta)};
      91         462 :     corners[i] = center + corner;
      92             :   }
      93          77 : }
      94             : 
      95             : void
      96          77 : SCMTriDuctMeshGenerator::ductCrossSec(std::vector<Point> & cross_sec,
      97             :                                       const std::vector<Point> & corners,
      98             :                                       unsigned int nrings,
      99             :                                       Real pitch) const
     100             : {
     101          77 :   cross_sec.clear();
     102             : 
     103          77 :   if (pitch <= 0.0)
     104           0 :     mooseError("The 'pitch' must be > 0.");
     105             : 
     106          77 :   if (nrings < 2)
     107           0 :     mooseError("The 'nrings' must be >= 2.");
     108             : 
     109             :   // Side length of the regular hexagon (distance between adjacent corners)
     110          77 :   const Real side_length = (corners[0] - corners[1]).norm();
     111             : 
     112             :   // start_offset = (hexagon_side - (nrings - 2)*pitch) / 2
     113          77 :   const Real start_offset = 0.5 * (side_length - (nrings - 2) * pitch);
     114             : 
     115             :   // If this is negative, the requested layout cannot fit.
     116          77 :   if (start_offset < 0.0)
     117           0 :     mooseError("SCMTriDuctMeshGenerator: computed start_offset is negative (",
     118             :                start_offset,
     119             :                "). Check 'nrings', 'pitch', and duct side length.");
     120             : 
     121         539 :   for (size_t i = 0; i < corners.size(); i++)
     122             :   {
     123         462 :     const Point left = corners[i];
     124         462 :     const Point right = corners[(i + 1) % corners.size()];
     125             : 
     126             :     // Always include the left corner
     127         462 :     cross_sec.push_back(left);
     128             : 
     129             :     // Unit direction along this side
     130         462 :     const Point direc = (right - left).unit();
     131             : 
     132             :     // Add exactly (nrings - 1) points after the left corner
     133             :     // at start_offset + k*pitch, k = 0..(nrings-2)
     134        2220 :     for (const auto k : make_range(nrings - 1))
     135             :     {
     136        1758 :       const Real offset_from_corner = start_offset + k * pitch;
     137        1758 :       cross_sec.push_back(left + direc * offset_from_corner);
     138             :     }
     139             :   }
     140          77 : }
     141             : 
     142             : size_t
     143      197970 : SCMTriDuctMeshGenerator::ductPointIndex(unsigned int points_per_layer,
     144             :                                         unsigned int layer,
     145             :                                         unsigned int point) const
     146             : {
     147      197970 :   return layer * points_per_layer + point;
     148             : }
     149             : 
     150             : void
     151          77 : SCMTriDuctMeshGenerator::ductPoints(std::vector<Point> & points,
     152             :                                     const std::vector<Point> & cross_sec,
     153             :                                     const std::vector<Real> & z_layers) const
     154             : {
     155          77 :   points.resize(cross_sec.size() * z_layers.size());
     156        1339 :   for (size_t i = 0; i < z_layers.size(); i++)
     157       42632 :     for (size_t j = 0; j < cross_sec.size(); j++)
     158       41370 :       points[ductPointIndex(cross_sec.size(), i, j)] =
     159       41370 :           Point(cross_sec[j](0), cross_sec[j](1), z_layers[i]);
     160          77 : }
     161             : 
     162             : void
     163          77 : SCMTriDuctMeshGenerator::ductElems(std::vector<std::vector<size_t>> & elem_point_indices,
     164             :                                    unsigned int n_layers,
     165             :                                    unsigned int points_per_layer) const
     166             : {
     167          77 :   elem_point_indices.clear();
     168        1262 :   for (unsigned int i = 0; i < n_layers - 1; i++)
     169             :   {
     170             :     unsigned int bottom = i;
     171        1185 :     unsigned int top = i + 1;
     172       40335 :     for (unsigned int j = 0; j < points_per_layer; j++)
     173             :     {
     174             :       unsigned int left = j;
     175       39150 :       unsigned int right = (j + 1) % points_per_layer;
     176      117450 :       elem_point_indices.push_back({ductPointIndex(points_per_layer, bottom, left),
     177       39150 :                                     ductPointIndex(points_per_layer, bottom, right),
     178       39150 :                                     ductPointIndex(points_per_layer, top, right),
     179       39150 :                                     ductPointIndex(points_per_layer, top, left)});
     180             :     }
     181             :   }
     182          77 : }
     183             : 
     184             : void
     185          77 : SCMTriDuctMeshGenerator::buildDuct(std::unique_ptr<MeshBase> & mesh,
     186             :                                    std::vector<Node *> & duct_nodes,
     187             :                                    const std::vector<Point> & points,
     188             :                                    const std::vector<std::vector<size_t>> & elem_point_indices,
     189             :                                    SubdomainID block) const
     190             : {
     191             :   // Create mesh nodes for all duct points and keep a local index -> Node* map
     192          77 :   duct_nodes.clear();
     193          77 :   duct_nodes.reserve(points.size());
     194       41447 :   for (const auto & p : points)
     195       41370 :     duct_nodes.push_back(mesh->add_point(p));
     196             : 
     197             :   // Create QUAD4 surface elements using libMesh factory style
     198       39227 :   for (const auto & elem_indices : elem_point_indices)
     199             :   {
     200             :     mooseAssert(elem_indices.size() == 4,
     201             :                 "Expected 4 node indices per element when building QUAD4 elements.");
     202             : 
     203       39150 :     auto elem = Elem::build(ElemType::QUAD4);
     204       39150 :     elem->subdomain_id() = block;
     205             : 
     206             :     // Set the 4 nodes of the QUAD4
     207      195750 :     for (unsigned int i = 0; i < 4; ++i)
     208             :     {
     209      156600 :       const auto idx = elem_indices[i];
     210             :       mooseAssert(idx < duct_nodes.size(), "Element node index out of range.");
     211      156600 :       elem->set_node(i, duct_nodes[idx]);
     212             :     }
     213             : 
     214             :     // Hand ownership to the mesh
     215       39150 :     mesh->add_elem(elem.release());
     216       39150 :   }
     217          77 : }