13 #include "libmesh/cell_prism6.h" 14 #include "libmesh/unstructured_mesh.h" 18 DetailedQuadSubChannelMeshGenerator,
27 "Creates a detailed mesh of subchannels in a square lattice arrangement");
30 params.
addParam<
Real>(
"unheated_length_entry", 0.0,
"Unheated length at entry [m]");
32 params.
addParam<
Real>(
"unheated_length_exit", 0.0,
"Unheated length at exit [m]");
33 params.
addRequiredParam<
unsigned int>(
"n_cells",
"The number of cells in the axial direction");
34 params.
addRequiredParam<
unsigned int>(
"nx",
"Number of channels in the x direction [-]");
35 params.
addRequiredParam<
unsigned int>(
"ny",
"Number of channels in the y direction [-]");
38 "The side gap, not to be confused with the gap between pins, this refers to the gap " 39 "next to the duct or else the distance between the subchannel centroid to the duct wall." 40 "distance(edge pin center, duct wall) = pitch / 2 + side_gap [m]");
42 params.
addParam<
unsigned int>(
"block_id", 0,
"Block ID used for the mesh subdomain.");
49 _unheated_length_entry(getParam<
Real>(
"unheated_length_entry")),
50 _heated_length(getParam<
Real>(
"heated_length")),
51 _unheated_length_exit(getParam<
Real>(
"unheated_length_exit")),
52 _pitch(getParam<
Real>(
"pitch")),
53 _pin_diameter(getParam<
Real>(
"pin_diameter")),
54 _n_cells(getParam<unsigned
int>(
"n_cells")),
55 _nx(getParam<unsigned
int>(
"nx")),
56 _ny(getParam<unsigned
int>(
"ny")),
57 _n_channels(_nx * _ny),
58 _side_gap(getParam<
Real>(
"side_gap")),
59 _block_id(getParam<unsigned
int>(
"block_id"))
63 for (
unsigned int i = 0; i <
_n_cells + 1; i++)
70 for (
unsigned int j = 0;
j < 3;
j++)
77 for (
unsigned int iy = 0; iy <
_ny; iy++)
79 for (
unsigned int ix = 0; ix <
_nx; ix++)
81 unsigned int i_ch =
_nx * iy + ix;
82 bool is_corner = (ix == 0 && iy == 0) || (ix ==
_nx - 1 && iy == 0) ||
83 (ix == 0 && iy ==
_ny - 1) || (ix ==
_nx - 1 && iy ==
_ny - 1);
84 bool is_edge = (ix == 0 || iy == 0 || ix ==
_nx - 1 || iy ==
_ny - 1);
102 std::unique_ptr<MeshBase>
106 BoundaryInfo & boundary_info = mesh_base->get_boundary_info();
107 mesh_base->set_spatial_dimension(3);
110 const unsigned int theta_res = 16;
112 const unsigned int points_per_quad = theta_res / 4 + 1;
117 const unsigned int points_per_center = points_per_quad * 4 + 1;
120 const unsigned int points_per_corner = points_per_quad * 1 + 1 + 3;
123 const unsigned int points_per_side = points_per_quad * 2 + 1 + 2;
127 const unsigned int elems_per_center = theta_res + 4;
128 const unsigned int elems_per_corner = theta_res / 4 + 4;
129 const unsigned int elems_per_side = theta_res / 2 + 4;
132 unsigned int n_center, n_side, n_corner;
148 n_side = 2 * (
_ny - 2) + 2 * (
_nx - 2);
153 const unsigned int points_per_level =
154 n_corner * points_per_corner + n_side * points_per_side + n_center * points_per_center;
155 const unsigned int elems_per_level =
156 n_corner * elems_per_corner + n_side * elems_per_side + n_center * elems_per_center;
157 const unsigned int n_points = points_per_level * (
_n_cells + 1);
158 const unsigned int n_elems = elems_per_level *
_n_cells;
159 mesh_base->reserve_nodes(n_points);
160 mesh_base->reserve_elem(n_elems);
163 std::array<Point, theta_res + 1> circle_points;
166 for (
unsigned int i = 0; i < theta_res + 1; i++)
168 circle_points[i](0) =
radius * std::cos(theta);
169 circle_points[i](1) =
radius * std::sin(theta);
170 theta += 2 * M_PI / theta_res;
178 const Real shrink_factor = 0.99999;
179 std::array<Point, 4> quadrant_centers;
180 quadrant_centers[0] = Point(
_pitch * 0.5 * shrink_factor,
_pitch * 0.5 * shrink_factor, 0);
181 quadrant_centers[1] = Point(-
_pitch * 0.5 * shrink_factor,
_pitch * 0.5 * shrink_factor, 0);
182 quadrant_centers[2] = Point(-
_pitch * 0.5 * shrink_factor, -
_pitch * 0.5 * shrink_factor, 0);
183 quadrant_centers[3] = Point(
_pitch * 0.5 * shrink_factor, -
_pitch * 0.5 * shrink_factor, 0);
185 const unsigned int m = theta_res / 4;
193 std::array<Point, points_per_center> center_points;
196 for (
unsigned int i = 0; i < 4; i++)
206 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
208 auto c_pt = circle_points[start - ii];
209 center_points[i * points_per_quad + ii + 1] = quadrant_centers[i] + c_pt;
221 std::array<Point, points_per_corner> tl_corner_points;
223 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
225 auto c_pt = circle_points[2 * m - ii];
226 tl_corner_points[ii + 1] = quadrant_centers[3] + c_pt;
228 tl_corner_points[points_per_quad + 1] = Point(
_pitch * 0.5 * shrink_factor,
_side_gap, 0);
230 tl_corner_points[points_per_quad + 3] = Point(-
_side_gap, -
_pitch * 0.5 * shrink_factor, 0);
240 std::array<Point, points_per_corner> tr_corner_points;
242 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
244 auto c_pt = circle_points[m - ii];
245 tr_corner_points[ii + 1] = quadrant_centers[2] + c_pt;
247 tr_corner_points[points_per_quad + 1] = Point(
_side_gap, -
_pitch * 0.5 * shrink_factor, 0);
249 tr_corner_points[points_per_quad + 3] = Point(-
_pitch * 0.5 * shrink_factor,
_side_gap, 0);
259 std::array<Point, points_per_corner> bl_corner_points;
261 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
263 auto c_pt = circle_points[3 * m - ii];
264 bl_corner_points[ii + 1] = quadrant_centers[0] + c_pt;
266 bl_corner_points[points_per_quad + 1] = Point(-
_side_gap,
_pitch * 0.5 * shrink_factor, 0);
268 bl_corner_points[points_per_quad + 3] = Point(
_pitch * 0.5 * shrink_factor, -
_side_gap, 0);
278 std::array<Point, points_per_corner> br_corner_points;
280 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
282 auto c_pt = circle_points[4 * m - ii];
283 br_corner_points[ii + 1] = quadrant_centers[1] + c_pt;
285 br_corner_points[points_per_quad + 1] = Point(-
_pitch * 0.5 * shrink_factor, -
_side_gap, 0);
287 br_corner_points[points_per_quad + 3] = Point(
_side_gap,
_pitch * 0.5 * shrink_factor, 0);
297 std::array<Point, points_per_side> top_points;
299 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
301 auto c_pt = circle_points[m - ii];
302 top_points[ii + 1] = quadrant_centers[2] + c_pt;
304 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
306 auto c_pt = circle_points[2 * m - ii];
307 top_points[points_per_quad + ii + 1] = quadrant_centers[3] + c_pt;
309 top_points[2 * points_per_quad + 1] = Point(
_pitch * 0.5 * shrink_factor,
_side_gap, 0);
310 top_points[2 * points_per_quad + 2] = Point(-
_pitch * 0.5 * shrink_factor,
_side_gap, 0);
320 std::array<Point, points_per_side> left_points;
322 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
324 auto c_pt = circle_points[2 * m - ii];
325 left_points[ii + 1] = quadrant_centers[3] + c_pt;
327 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
329 auto c_pt = circle_points[3 * m - ii];
330 left_points[points_per_quad + ii + 1] = quadrant_centers[0] + c_pt;
332 left_points[2 * points_per_quad + 1] = Point(-
_side_gap,
_pitch * 0.5 * shrink_factor, 0);
333 left_points[2 * points_per_quad + 2] = Point(-
_side_gap, -
_pitch * 0.5 * shrink_factor, 0);
343 std::array<Point, points_per_side> bottom_points;
345 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
347 auto c_pt = circle_points[3 * m - ii];
348 bottom_points[ii + 1] = quadrant_centers[0] + c_pt;
350 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
352 auto c_pt = circle_points[4 * m - ii];
353 bottom_points[points_per_quad + ii + 1] = quadrant_centers[1] + c_pt;
355 bottom_points[2 * points_per_quad + 1] = Point(-
_pitch * 0.5 * shrink_factor, -
_side_gap, 0);
356 bottom_points[2 * points_per_quad + 2] = Point(
_pitch * 0.5 * shrink_factor, -
_side_gap, 0);
366 std::array<Point, points_per_side> right_points;
368 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
370 auto c_pt = circle_points[4 * m - ii];
371 right_points[ii + 1] = quadrant_centers[1] + c_pt;
373 for (
unsigned int ii = 0; ii < points_per_quad; ii++)
375 auto c_pt = circle_points[m - ii];
376 right_points[points_per_quad + ii + 1] = quadrant_centers[2] + c_pt;
378 right_points[2 * points_per_quad + 1] = Point(
_side_gap, -
_pitch * 0.5 * shrink_factor, 0);
379 right_points[2 * points_per_quad + 2] = Point(
_side_gap,
_pitch * 0.5 * shrink_factor, 0);
385 unsigned int node_id = 0;
388 for (
unsigned int iy = 0; iy <
_ny; iy++)
390 Point y0 = {0,
_pitch * iy - offset_y, 0};
391 for (
unsigned int ix = 0; ix <
_nx; ix++)
393 Point x0 = {
_pitch * ix - offset_x, 0, 0};
397 if (
_nx == 1 && iy == 0)
399 for (
unsigned int i = 0; i < points_per_side; i++)
400 mesh_base->add_point(bottom_points[i] + x0 + y0 + z0, node_id++);
402 if (
_nx == 1 && iy == 1)
404 for (
unsigned int i = 0; i < points_per_side; i++)
405 mesh_base->add_point(top_points[i] + x0 + y0 + z0, node_id++);
407 if (
_ny == 1 && ix == 0)
409 for (
unsigned int i = 0; i < points_per_side; i++)
410 mesh_base->add_point(left_points[i] + x0 + y0 + z0, node_id++);
412 if (
_ny == 1 && ix == 1)
414 for (
unsigned int i = 0; i < points_per_side; i++)
415 mesh_base->add_point(right_points[i] + x0 + y0 + z0, node_id++);
423 unsigned int node_id = 0;
426 for (
unsigned int iy = 0; iy <
_ny; iy++)
428 Point y0 = {0,
_pitch * iy - offset_y, 0};
429 for (
unsigned int ix = 0; ix <
_nx; ix++)
431 Point x0 = {
_pitch * ix - offset_x, 0, 0};
439 for (
unsigned int i = 0; i < points_per_side; i++)
440 mesh_base->add_point(bottom_points[i] + x0 + y0 + z0, node_id++);
442 else if (iy ==
_ny - 1)
444 for (
unsigned int i = 0; i < points_per_side; i++)
445 mesh_base->add_point(top_points[i] + x0 + y0 + z0, node_id++);
449 for (
unsigned int i = 0; i < points_per_center; i++)
450 mesh_base->add_point(center_points[i] + x0 + y0 + z0, node_id++);
457 for (
unsigned int i = 0; i < points_per_side; i++)
458 mesh_base->add_point(left_points[i] + x0 + y0 + z0, node_id++);
460 else if (ix ==
_nx - 1)
462 for (
unsigned int i = 0; i < points_per_side; i++)
463 mesh_base->add_point(right_points[i] + x0 + y0 + z0, node_id++);
467 for (
unsigned int i = 0; i < points_per_center; i++)
468 mesh_base->add_point(center_points[i] + x0 + y0 + z0, node_id++);
477 unsigned int node_id = 0;
480 for (
unsigned int iy = 0; iy <
_ny; iy++)
482 Point y0 = {0,
_pitch * iy - offset_y, 0};
483 for (
unsigned int ix = 0; ix <
_nx; ix++)
485 Point x0 = {
_pitch * ix - offset_x, 0, 0};
486 if (ix == 0 && iy == 0)
491 for (
unsigned int i = 0; i < points_per_corner; i++)
492 mesh_base->add_point(bl_corner_points[i] + x0 + y0 + z0, node_id++);
495 else if (ix ==
_nx - 1 && iy == 0)
500 for (
unsigned int i = 0; i < points_per_corner; i++)
501 mesh_base->add_point(br_corner_points[i] + x0 + y0 + z0, node_id++);
504 else if (ix == 0 && iy ==
_ny - 1)
509 for (
unsigned int i = 0; i < points_per_corner; i++)
510 mesh_base->add_point(tl_corner_points[i] + x0 + y0 + z0, node_id++);
513 else if (ix ==
_nx - 1 && iy ==
_ny - 1)
518 for (
unsigned int i = 0; i < points_per_corner; i++)
519 mesh_base->add_point(tr_corner_points[i] + x0 + y0 + z0, node_id++);
522 else if (ix == 0 && (iy !=
_ny - 1 || iy != 0))
527 for (
unsigned int i = 0; i < points_per_side; i++)
528 mesh_base->add_point(left_points[i] + x0 + y0 + z0, node_id++);
531 else if (ix ==
_nx - 1 && (iy !=
_ny - 1 || iy != 0))
536 for (
unsigned int i = 0; i < points_per_side; i++)
537 mesh_base->add_point(right_points[i] + x0 + y0 + z0, node_id++);
540 else if (iy == 0 && (ix !=
_nx - 1 || ix != 0))
545 for (
unsigned int i = 0; i < points_per_side; i++)
546 mesh_base->add_point(bottom_points[i] + x0 + y0 + z0, node_id++);
549 else if (iy ==
_ny - 1 && (ix !=
_nx - 1 || ix != 0))
554 for (
unsigned int i = 0; i < points_per_side; i++)
555 mesh_base->add_point(top_points[i] + x0 + y0 + z0, node_id++);
563 for (
unsigned int i = 0; i < points_per_center; i++)
564 mesh_base->add_point(center_points[i] + x0 + y0 + z0, node_id++);
576 unsigned int elem_id = 0;
577 for (
unsigned int iy = 0; iy <
_ny; iy++)
579 for (
unsigned int ix = 0; ix <
_nx; ix++)
581 unsigned int i_ch =
_nx * iy + ix;
582 for (
unsigned int iz = 0; iz <
_n_cells; iz++)
584 for (
unsigned int i = 0; i < elems_per_side; i++)
586 Elem * elem =
new Prism6;
588 elem->set_id(elem_id++);
589 elem = mesh_base->add_elem(elem);
591 unsigned int indx1 = iz * points_per_side + points_per_side * (
_n_cells + 1) * i_ch;
594 (iz + 1) * points_per_side + points_per_side * (
_n_cells + 1) * i_ch;
595 unsigned int elems_per_channel = elems_per_side;
596 elem->set_node(0, mesh_base->node_ptr(indx1));
597 elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
598 if (i != elems_per_channel - 1)
599 elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
601 elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
603 elem->set_node(3, mesh_base->node_ptr(indx2));
604 elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
605 if (i != elems_per_channel - 1)
606 elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
608 elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
611 boundary_info.add_side(elem, 0, 0);
613 boundary_info.add_side(elem, 4, 1);
618 boundary_info.sideset_name(0) =
"inlet";
619 boundary_info.sideset_name(1) =
"outlet";
621 mesh_base->prepare_for_use();
625 unsigned int elem_id = 0;
626 unsigned int number_of_corner = 0;
627 unsigned int number_of_side = 0;
628 unsigned int number_of_center = 0;
629 unsigned int elems_per_channel = 0;
630 unsigned int points_per_channel = 0;
631 for (
unsigned int iy = 0; iy <
_ny; iy++)
633 for (
unsigned int ix = 0; ix <
_nx; ix++)
635 unsigned int i_ch =
_nx * iy + ix;
641 elems_per_channel = elems_per_side;
642 points_per_channel = points_per_side;
648 elems_per_channel = elems_per_center;
649 points_per_channel = points_per_center;
651 for (
unsigned int iz = 0; iz <
_n_cells; iz++)
653 unsigned int elapsed_points = number_of_corner * points_per_corner * (
_n_cells + 1) +
654 number_of_side * points_per_side * (
_n_cells + 1) +
655 number_of_center * points_per_center * (
_n_cells + 1) -
656 points_per_channel * (
_n_cells + 1);
658 unsigned int indx1 = iz * points_per_channel + elapsed_points;
660 unsigned int indx2 = (iz + 1) * points_per_channel + elapsed_points;
662 for (
unsigned int i = 0; i < elems_per_channel; i++)
664 Elem * elem =
new Prism6;
666 elem->set_id(elem_id++);
667 elem = mesh_base->add_elem(elem);
669 elem->set_node(0, mesh_base->node_ptr(indx1));
670 elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
671 if (i != elems_per_channel - 1)
672 elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
674 elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
676 elem->set_node(3, mesh_base->node_ptr(indx2));
677 elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
678 if (i != elems_per_channel - 1)
679 elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
681 elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
684 boundary_info.add_side(elem, 0, 0);
686 boundary_info.add_side(elem, 4, 1);
691 boundary_info.sideset_name(0) =
"inlet";
692 boundary_info.sideset_name(1) =
"outlet";
694 mesh_base->prepare_for_use();
698 unsigned int elem_id = 0;
699 unsigned int number_of_corner = 0;
700 unsigned int number_of_side = 0;
701 unsigned int number_of_center = 0;
702 unsigned int elems_per_channel = 0;
703 unsigned int points_per_channel = 0;
704 for (
unsigned int iy = 0; iy <
_ny; iy++)
706 for (
unsigned int ix = 0; ix <
_nx; ix++)
708 unsigned int i_ch =
_nx * iy + ix;
713 elems_per_channel = elems_per_corner;
714 points_per_channel = points_per_corner;
719 elems_per_channel = elems_per_side;
720 points_per_channel = points_per_side;
725 elems_per_channel = elems_per_center;
726 points_per_channel = points_per_center;
728 for (
unsigned int iz = 0; iz <
_n_cells; iz++)
730 unsigned int elapsed_points = number_of_corner * points_per_corner * (
_n_cells + 1) +
731 number_of_side * points_per_side * (
_n_cells + 1) +
732 number_of_center * points_per_center * (
_n_cells + 1) -
733 points_per_channel * (
_n_cells + 1);
735 unsigned int indx1 = iz * points_per_channel + elapsed_points;
737 unsigned int indx2 = (iz + 1) * points_per_channel + elapsed_points;
739 for (
unsigned int i = 0; i < elems_per_channel; i++)
741 Elem * elem =
new Prism6;
743 elem->set_id(elem_id++);
744 elem = mesh_base->add_elem(elem);
746 elem->set_node(0, mesh_base->node_ptr(indx1));
747 elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
748 if (i != elems_per_channel - 1)
749 elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
751 elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
753 elem->set_node(3, mesh_base->node_ptr(indx2));
754 elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
755 if (i != elems_per_channel - 1)
756 elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
758 elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
761 boundary_info.add_side(elem, 0, 0);
763 boundary_info.add_side(elem, 4, 1);
768 boundary_info.sideset_name(0) =
"inlet";
769 boundary_info.sideset_name(1) =
"outlet";
771 mesh_base->prepare_for_use();
registerMooseObject("SubChannelApp", SCMDetailedQuadSubChannelMeshGenerator)
registerMooseObjectRenamed("SubChannelApp", DetailedQuadSubChannelMeshGenerator, "06/30/2025 24:00", SCMDetailedQuadSubChannelMeshGenerator)
const Real _pitch
Distance between the neighbor fuel pins, pitch.
const unsigned int _n_cells
Number of cells in the axial direction.
const Real _heated_length
heated length of the fuel Pin
std::vector< Real > _z_grid
axial location of nodes
unsigned int _n_channels
Total number of subchannels.
const unsigned int _nx
Number of subchannels in the x direction.
const std::string & name() const
std::vector< std::vector< Real > > _subchannel_position
x,y coordinates of the subchannel centroids
const Real _side_gap
The side gap, not to be confused with the gap between pins, this refers to the gap next to the duct o...
static InputParameters validParams()
const Real _pin_diameter
fuel Pin diameter
SCMDetailedQuadSubChannelMeshGenerator(const InputParameters ¶meters)
const unsigned int _ny
Number of subchannels in the y direction.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Mesh generator that builds a 3D mesh representing quadrilateral subchannels.
const Real _unheated_length_exit
unheated length of the fuel Pin at the exit of the assembly
std::vector< EChannelType > _subch_type
Subchannel type.
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
EChannelType getSubchannelType(unsigned int index) const
std::unique_ptr< MeshBase > buildMeshBaseObject(unsigned int dim=libMesh::invalid_uint)
const Real _unheated_length_entry
unheated length of the fuel Pin at the entry of the assembly
const unsigned int & _block_id
Subdomain ID used for the mesh block.
static InputParameters validParams()
virtual std::unique_ptr< MeshBase > generate() override
void ErrorVector unsigned int