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 "SCMDetailedTriAssemblyMeshGenerator.h"
11 : #include "TriSubChannelMesh.h"
12 : #include <array>
13 : #include <cmath>
14 : #include <memory>
15 : #include "libmesh/cell_prism6.h"
16 : #include "libmesh/unstructured_mesh.h"
17 :
18 : registerMooseObject("SubChannelApp", SCMDetailedTriAssemblyMeshGenerator);
19 : registerMooseObjectRenamed("SubChannelApp",
20 : SCMDetailedTriSubChannelMeshGenerator,
21 : "06/30/2027 24:00",
22 : SCMDetailedTriAssemblyMeshGenerator);
23 : registerMooseObjectRenamed("SubChannelApp",
24 : DetailedTriSubChannelMeshGenerator,
25 : "06/30/2027 24:00",
26 : SCMDetailedTriAssemblyMeshGenerator);
27 : registerMooseObjectRenamed("SubChannelApp",
28 : SCMDetailedTriPinMeshGenerator,
29 : "06/30/2027 24:00",
30 : SCMDetailedTriAssemblyMeshGenerator);
31 : registerMooseObjectRenamed("SubChannelApp",
32 : DetailedTriPinMeshGenerator,
33 : "06/30/2027 24:00",
34 : SCMDetailedTriAssemblyMeshGenerator);
35 :
36 : InputParameters
37 66 : SCMDetailedTriAssemblyMeshGenerator::validParams()
38 : {
39 66 : InputParameters params = MeshGenerator::validParams();
40 66 : params.addClassDescription(
41 : "Creates a detailed mesh of subchannels and pins in a triangular lattice arrangement");
42 132 : params.addRequiredParam<Real>("pitch", "Pitch [m]");
43 132 : params.addRequiredParam<Real>("pin_diameter", "Rod diameter [m]");
44 132 : params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
45 132 : params.addRequiredParam<Real>("heated_length", "Heated length [m]");
46 132 : params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");
47 132 : params.addRequiredParam<unsigned int>(
48 : "nrings",
49 : "Number of fuel-pin rings per assembly, counting the center pin as the first ring [-]");
50 132 : params.addRequiredParam<Real>("flat_to_flat",
51 : "Flat to flat distance for the hexagonal assembly [m]");
52 132 : params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
53 198 : params.addRangeCheckedParam<unsigned int>("num_sectors",
54 132 : 16,
55 : "num_sectors>=4",
56 : "Number of azimuthal sectors used to discretize each "
57 : "circular pin cross section.");
58 132 : params.addParam<unsigned int>("block_id", 0, "Subchannel block id.");
59 132 : params.deprecateParam("block_id", "subchannel_block_id", "07/01/2027");
60 132 : params.addParam<unsigned int>("pin_block_id", 1, "Fuel pin block id.");
61 132 : params.addParam<bool>("verbose_flag", false, "Flag to print out the mesh coordinates.");
62 66 : return params;
63 0 : }
64 :
65 33 : SCMDetailedTriAssemblyMeshGenerator::SCMDetailedTriAssemblyMeshGenerator(
66 33 : const InputParameters & parameters)
67 : : MeshGenerator(parameters),
68 33 : _unheated_length_entry(getParam<Real>("unheated_length_entry")),
69 66 : _heated_length(getParam<Real>("heated_length")),
70 66 : _unheated_length_exit(getParam<Real>("unheated_length_exit")),
71 66 : _pitch(getParam<Real>("pitch")),
72 66 : _pin_diameter(getParam<Real>("pin_diameter")),
73 66 : _n_rings(getParam<unsigned int>("nrings")),
74 66 : _flat_to_flat(getParam<Real>("flat_to_flat")),
75 66 : _num_sectors(getParam<unsigned int>("num_sectors")),
76 66 : _subchannel_block_id(getParam<unsigned int>("subchannel_block_id")),
77 66 : _pin_block_id(getParam<unsigned int>("pin_block_id")),
78 66 : _n_cells(getParam<unsigned int>("n_cells")),
79 33 : _nrods(0),
80 33 : _n_channels(0),
81 66 : _verbose(getParam<bool>("verbose_flag")),
82 33 : _elem_id(0)
83 : {
84 33 : const Real L = _unheated_length_entry + _heated_length + _unheated_length_exit;
85 :
86 33 : if (_n_rings < 2)
87 0 : paramError("nrings",
88 : "'nrings' must be at least 2. In this mesh generator, the center pin counts as "
89 : "the first ring, so a 7-pin bundle uses nrings = 2.");
90 :
91 33 : if (_n_cells == 0)
92 0 : paramError("n_cells", "The number of axial cells must be greater than zero");
93 :
94 33 : if (L <= 0.0)
95 0 : mooseError("Total bundle length must be greater than zero");
96 :
97 33 : Real dz = L / _n_cells;
98 1266 : for (unsigned int i = 0; i < _n_cells + 1; i++)
99 1233 : _z_grid.push_back(dz * i);
100 :
101 : // x coordinate for the first position
102 : Real x0 = 0.0;
103 : // y coordinate for the first position
104 : Real y0 = 0.0;
105 : // x coordinate for the second position
106 : Real x1 = 0.0;
107 : // y coordinate for the second position dummy variable
108 : Real y1 = 0.0;
109 : // dummy variable
110 : Real a1 = 0.0;
111 : // dummy variable
112 : Real a2 = 0.0;
113 : // average x coordinate
114 : Real avg_coor_x = 0.0;
115 : // average y coordinate
116 : Real avg_coor_y = 0.0;
117 : // distance between two points
118 : Real dist = 0.0;
119 : // distance between two points
120 : Real dist0 = 0.0;
121 : // the indicator used while setting _gap_to_chan_map array
122 : std::vector<std::pair<unsigned int, unsigned int>> gap_fill;
123 33 : TriSubChannelMesh::pinPositions(_pin_position, _n_rings, _pitch, Point(0, 0));
124 33 : _nrods = _pin_position.size();
125 : // assign the pins to the corresponding rings
126 : unsigned int k = 0; // initialize the fuel Pin counter index
127 33 : _pins_in_rings.resize(_n_rings);
128 33 : _pins_in_rings[0].push_back(k++);
129 105 : for (unsigned int i = 1; i < _n_rings; i++)
130 774 : for (unsigned int j = 0; j < i * 6; j++)
131 702 : _pins_in_rings[i].push_back(k++);
132 : // Given the number of pins and number of fuel Pin rings, the number of subchannels can be
133 : // computed as follows:
134 : unsigned int chancount = 0;
135 : // Summing internal channels
136 105 : for (unsigned int j = 0; j < _n_rings - 1; j++)
137 72 : chancount += j * 6;
138 : // Adding external channels to the total count
139 33 : _n_channels = chancount + _nrods - 1 + (_n_rings - 1) * 6 + 6;
140 :
141 : // Utils for building the mesh
142 33 : _chan_to_pin_map.resize(_n_channels);
143 33 : _subch_type.resize(_n_channels);
144 33 : _subchannel_position.resize(_n_channels);
145 :
146 1635 : for (unsigned int i = 0; i < _n_channels; i++)
147 : {
148 1602 : _subchannel_position[i].reserve(3);
149 6408 : for (unsigned int j = 0; j < 3; j++)
150 : {
151 4806 : _subchannel_position.at(i).push_back(0.0);
152 : }
153 : }
154 :
155 : // create the subchannels
156 : k = 0; // initialize the subchannel counter index
157 105 : for (unsigned int i = 1; i < _n_rings; i++)
158 : {
159 : // find the closest Pin at back ring
160 774 : for (unsigned int j = 0; j < _pins_in_rings[i].size(); j++)
161 : {
162 702 : if (j == _pins_in_rings[i].size() - 1)
163 : {
164 72 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
165 72 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][0]);
166 72 : avg_coor_x =
167 72 : 0.5 * (_pin_position[_pins_in_rings[i][j]](0) + _pin_position[_pins_in_rings[i][0]](0));
168 72 : avg_coor_y =
169 72 : 0.5 * (_pin_position[_pins_in_rings[i][j]](1) + _pin_position[_pins_in_rings[i][0]](1));
170 : }
171 : else
172 : {
173 630 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
174 630 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j + 1]);
175 630 : avg_coor_x = 0.5 * (_pin_position[_pins_in_rings[i][j]](0) +
176 630 : _pin_position[_pins_in_rings[i][j + 1]](0));
177 630 : avg_coor_y = 0.5 * (_pin_position[_pins_in_rings[i][j]](1) +
178 630 : _pin_position[_pins_in_rings[i][j + 1]](1));
179 : }
180 : dist0 = 1.0e+5;
181 702 : _chan_to_pin_map[k].push_back(_pins_in_rings[i - 1][0]);
182 4572 : for (unsigned int l = 0; l < _pins_in_rings[i - 1].size(); l++)
183 : {
184 3870 : dist = std::sqrt(pow(_pin_position[_pins_in_rings[i - 1][l]](0) - avg_coor_x, 2) +
185 3870 : pow(_pin_position[_pins_in_rings[i - 1][l]](1) - avg_coor_y, 2));
186 :
187 3870 : if (dist < dist0)
188 : {
189 1710 : _chan_to_pin_map[k][2] = _pins_in_rings[i - 1][l];
190 : dist0 = dist;
191 : }
192 : }
193 702 : _subch_type[k] = EChannelType::CENTER;
194 702 : _orientation_map.insert(std::make_pair(k, 0.0));
195 702 : k = k + 1;
196 : }
197 :
198 : // find the closest Pin at front ring
199 774 : for (unsigned int j = 0; j < _pins_in_rings[i].size(); j++)
200 : {
201 702 : if (j == _pins_in_rings[i].size() - 1)
202 : {
203 72 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
204 72 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][0]);
205 72 : avg_coor_x =
206 72 : 0.5 * (_pin_position[_pins_in_rings[i][j]](0) + _pin_position[_pins_in_rings[i][0]](0));
207 72 : avg_coor_y =
208 72 : 0.5 * (_pin_position[_pins_in_rings[i][j]](1) + _pin_position[_pins_in_rings[i][0]](1));
209 : }
210 : else
211 : {
212 630 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
213 630 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j + 1]);
214 630 : avg_coor_x = 0.5 * (_pin_position[_pins_in_rings[i][j]](0) +
215 630 : _pin_position[_pins_in_rings[i][j + 1]](0));
216 630 : avg_coor_y = 0.5 * (_pin_position[_pins_in_rings[i][j]](1) +
217 630 : _pin_position[_pins_in_rings[i][j + 1]](1));
218 : }
219 : // if the outermost ring, set the edge subchannels first... then the corner subchannels
220 702 : if (i == _n_rings - 1)
221 : {
222 : // add edges
223 432 : _subch_type[k] = EChannelType::EDGE; // an edge subchannel is created
224 432 : k = k + 1;
225 432 : if (j % i == 0)
226 : {
227 : // corner subchannel
228 198 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
229 198 : _subch_type[k] = EChannelType::CORNER;
230 198 : k = k + 1;
231 : }
232 : // if not the outer most ring
233 : }
234 : else
235 : {
236 : dist0 = 1.0e+5;
237 270 : _chan_to_pin_map[k].push_back(_pins_in_rings[i + 1][0]);
238 3942 : for (unsigned int l = 0; l < _pins_in_rings[i + 1].size(); l++)
239 : {
240 3672 : dist = std::sqrt(pow(_pin_position[_pins_in_rings[i + 1][l]](0) - avg_coor_x, 2) +
241 3672 : pow(_pin_position[_pins_in_rings[i + 1][l]](1) - avg_coor_y, 2));
242 3672 : if (dist < dist0)
243 : {
244 1323 : _chan_to_pin_map[k][2] = _pins_in_rings[i + 1][l];
245 : dist0 = dist;
246 : }
247 : }
248 270 : _subch_type[k] = EChannelType::CENTER;
249 270 : _orientation_map.insert(std::make_pair(k, libMesh::pi));
250 270 : k = k + 1;
251 : }
252 : }
253 : }
254 :
255 1635 : for (auto & pin : _chan_to_pin_map)
256 : pin.shrink_to_fit();
257 :
258 : // set the subchannel positions
259 33 : Real _duct_to_pin_gap =
260 33 : 0.5 * (_flat_to_flat - (_n_rings - 1) * _pitch * std::sqrt(3.0) - _pin_diameter);
261 1635 : for (unsigned int i = 0; i < _n_channels; i++)
262 : {
263 1602 : if (_subch_type[i] == EChannelType::CENTER)
264 : {
265 972 : _subchannel_position[i][0] =
266 972 : (_pin_position[_chan_to_pin_map[i][0]](0) + _pin_position[_chan_to_pin_map[i][1]](0) +
267 972 : _pin_position[_chan_to_pin_map[i][2]](0)) /
268 : 3.0;
269 972 : _subchannel_position[i][1] =
270 972 : (_pin_position[_chan_to_pin_map[i][0]](1) + _pin_position[_chan_to_pin_map[i][1]](1) +
271 972 : _pin_position[_chan_to_pin_map[i][2]](1)) /
272 : 3.0;
273 : }
274 630 : else if (_subch_type[i] == EChannelType::EDGE)
275 : {
276 22464 : for (unsigned int j = 0; j < _n_channels; j++)
277 : {
278 22032 : if (_subch_type[j] == EChannelType::CENTER &&
279 13608 : ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][0] &&
280 432 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][1]) ||
281 13176 : (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][1] &&
282 432 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][0])))
283 : {
284 432 : x0 = _pin_position[_chan_to_pin_map[j][2]](0);
285 432 : y0 = _pin_position[_chan_to_pin_map[j][2]](1);
286 : }
287 21600 : else if (_subch_type[j] == EChannelType::CENTER &&
288 13176 : ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][0] &&
289 0 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][2]) ||
290 13176 : (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][2] &&
291 234 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][0])))
292 : {
293 0 : x0 = _pin_position[_chan_to_pin_map[j][1]](0);
294 0 : y0 = _pin_position[_chan_to_pin_map[j][1]](1);
295 : }
296 21600 : else if (_subch_type[j] == EChannelType::CENTER &&
297 13176 : ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][1] &&
298 432 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][2]) ||
299 13176 : (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][2] &&
300 234 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][1])))
301 : {
302 0 : x0 = _pin_position[_chan_to_pin_map[j][0]](0);
303 0 : y0 = _pin_position[_chan_to_pin_map[j][0]](1);
304 : }
305 22032 : x1 = 0.5 *
306 22032 : (_pin_position[_chan_to_pin_map[i][0]](0) + _pin_position[_chan_to_pin_map[i][1]](0));
307 22032 : y1 = 0.5 *
308 22032 : (_pin_position[_chan_to_pin_map[i][0]](1) + _pin_position[_chan_to_pin_map[i][1]](1));
309 22032 : a1 = _pin_diameter / 2.0 + _duct_to_pin_gap / 2.0;
310 22032 : a2 = std::sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) + a1;
311 22032 : _subchannel_position[i][0] = (a2 * x1 - a1 * x0) / (a2 - a1);
312 22032 : _subchannel_position[i][1] = (a2 * y1 - a1 * y0) / (a2 - a1);
313 : } // j
314 : }
315 198 : else if (_subch_type[i] == EChannelType::CORNER)
316 : {
317 198 : x0 = _pin_position[0](0);
318 198 : y0 = _pin_position[0](1);
319 198 : x1 = _pin_position[_chan_to_pin_map[i][0]](0);
320 198 : y1 = _pin_position[_chan_to_pin_map[i][0]](1);
321 198 : a1 = _pin_diameter / 2.0 + _duct_to_pin_gap / 2.0;
322 198 : a2 = std::sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) + a1;
323 198 : _subchannel_position[i][0] = (a2 * x1 - a1 * x0) / (a2 - a1);
324 198 : _subchannel_position[i][1] = (a2 * y1 - a1 * y0) / (a2 - a1);
325 : }
326 : }
327 33 : }
328 :
329 : void
330 735 : SCMDetailedTriAssemblyMeshGenerator::generatePin(std::unique_ptr<MeshBase> & mesh_base,
331 : const Point & center)
332 : {
333 735 : const Real dalpha = 360. / _num_sectors;
334 735 : const Real radius = _pin_diameter / 2.;
335 :
336 : // Add a center node and radial boundary nodes on each axial level so each pin is discretized into
337 : // triangular prism sectors.
338 : std::vector<std::vector<Node *>> nodes;
339 735 : nodes.resize(_n_cells + 1);
340 26430 : for (unsigned int k = 0; k < _n_cells + 1; k++)
341 : {
342 25695 : const Real elev = _z_grid[k];
343 25695 : nodes[k].push_back(mesh_base->add_point(Point(center(0), center(1), elev)));
344 : Real alpha = 0.;
345 436815 : for (unsigned int i = 0; i < _num_sectors; i++, alpha += dalpha)
346 : {
347 411120 : const Real dx = radius * std::cos(alpha * M_PI / 180.);
348 411120 : const Real dy = radius * std::sin(alpha * M_PI / 180.);
349 411120 : nodes[k].push_back(mesh_base->add_point(Point(center(0) + dx, center(1) + dy, elev)));
350 : }
351 : }
352 :
353 : // Add the pin volume elements, linking matching radial sectors between adjacent axial levels.
354 25695 : for (unsigned int k = 0; k < _n_cells; k++)
355 424320 : for (unsigned int i = 0; i < _num_sectors; i++)
356 : {
357 399360 : Elem * elem = mesh_base->add_elem(std::make_unique<Prism6>());
358 399360 : elem->subdomain_id() = _pin_block_id;
359 399360 : elem->set_id(_elem_id++);
360 : const unsigned int ctr_idx = 0;
361 399360 : const unsigned int idx1 = (i % _num_sectors) + 1;
362 399360 : const unsigned int idx2 = ((i + 1) % _num_sectors) + 1;
363 399360 : elem->set_node(0, nodes[k][ctr_idx]);
364 399360 : elem->set_node(1, nodes[k][idx1]);
365 399360 : elem->set_node(2, nodes[k][idx2]);
366 399360 : elem->set_node(3, nodes[k + 1][ctr_idx]);
367 399360 : elem->set_node(4, nodes[k + 1][idx1]);
368 399360 : elem->set_node(5, nodes[k + 1][idx2]);
369 : }
370 735 : }
371 :
372 : std::unique_ptr<MeshBase>
373 33 : SCMDetailedTriAssemblyMeshGenerator::generate()
374 : {
375 33 : auto mesh_base = buildMeshBaseObject();
376 : BoundaryInfo & boundary_info = mesh_base->get_boundary_info();
377 33 : mesh_base->set_spatial_dimension(3);
378 : // Define the resolution (the number of points used to represent a circle).
379 : // This must be divisible by 4.
380 : const unsigned int theta_res_triangle = 24; // TODO: parameterize
381 : const unsigned int theta_res_square = 24; // TODO: parameterize
382 : // Compute the number of points needed to represent one sixth and quadrant of a circle.
383 : const unsigned int points_per_sixth = theta_res_triangle / 6 + 1;
384 : const unsigned int points_per_quadrant = theta_res_square / 4 + 1;
385 :
386 : // Compute the points needed to represent one axial cross-flow of a subchannel.
387 : // For the center subchannel (sc) there is one center point plus the points from 3 side
388 : // circles.
389 33 : const unsigned int points_per_center = points_per_sixth * 3 + 1;
390 : // For the corner sc there is one center point plus the points from 1 side circle plus 3
391 : // corners
392 33 : const unsigned int points_per_corner = points_per_sixth * 1 + 1 + 3;
393 : // For the side sc there is one center point plus the points from 2 intersecting circles plus 2
394 : // corners
395 33 : const unsigned int points_per_side = points_per_quadrant * 2 + 1 + 2;
396 :
397 33 : if (_verbose)
398 : {
399 7 : _console << "Points per center: " << points_per_center << std::endl;
400 7 : _console << "Points per side: " << points_per_side << std::endl;
401 7 : _console << "Points per corner: " << points_per_corner << std::endl;
402 : }
403 :
404 : // Compute the number of elements (Prism6) which combined base creates the sub-channel
405 : // cross-section
406 33 : const unsigned int elems_per_center = theta_res_triangle * 3 / 6 + 3; // TODO: check
407 33 : const unsigned int elems_per_corner = theta_res_triangle / 6 + 4;
408 33 : const unsigned int elems_per_side = 2 * theta_res_square / 4 + 4;
409 33 : if (_verbose)
410 : {
411 7 : _console << "Elems per center: " << elems_per_center << std::endl;
412 7 : _console << "Elems per side: " << elems_per_side << std::endl;
413 7 : _console << "Elems per corner: " << elems_per_corner << std::endl;
414 : }
415 :
416 : // specify number and type of sub-channel
417 33 : unsigned int n_corner = 6;
418 33 : unsigned int n_side = (_n_rings - 1) * 6;
419 33 : unsigned int n_center = _n_channels - n_side - n_corner;
420 33 : if (_verbose)
421 : {
422 7 : _console << "Centers: " << n_center << std::endl;
423 7 : _console << "Sides: " << n_side << std::endl;
424 7 : _console << "Corners: " << n_corner << std::endl;
425 : }
426 :
427 : // Compute the total number of points and elements.
428 33 : const unsigned int points_per_level =
429 33 : n_corner * points_per_corner + n_side * points_per_side + n_center * points_per_center;
430 33 : const unsigned int elems_per_level =
431 33 : n_corner * elems_per_corner + n_side * elems_per_side + n_center * elems_per_center;
432 33 : if (_verbose)
433 : {
434 7 : _console << "Points per level: " << points_per_level << std::endl;
435 7 : _console << "Elements per level: " << elems_per_level << std::endl;
436 : }
437 33 : const unsigned int n_points = points_per_level * (_n_cells + 1);
438 33 : const unsigned int n_elems = elems_per_level * _n_cells;
439 33 : const unsigned int n_pins = _pin_position.size();
440 33 : const unsigned int pin_points = n_pins > 0 ? (_n_cells + 1) * (_num_sectors + 1) * n_pins : 0;
441 33 : const unsigned int pin_elems = n_pins > 0 ? _n_cells * _num_sectors * n_pins : 0;
442 33 : if (_verbose)
443 : {
444 7 : _console << "Number of points: " << n_points << std::endl;
445 7 : _console << "Number of elements: " << n_elems << std::endl;
446 : }
447 33 : mesh_base->reserve_nodes(n_points + pin_points);
448 33 : mesh_base->reserve_elem(n_elems + pin_elems);
449 : // Build an array of points arranged in a circle on the xy-plane. (last and first node overlap)
450 : // We build for both the square discretization in the edges and the triangular discretization
451 : // within the mesh
452 33 : const Real radius = _pin_diameter / 2.0;
453 : std::array<Point, theta_res_square + 1> circle_points_square;
454 : {
455 : Real theta = 0;
456 858 : for (unsigned int i = 0; i < theta_res_square + 1; i++)
457 : {
458 825 : circle_points_square[i](0) = radius * std::cos(theta);
459 825 : circle_points_square[i](1) = radius * std::sin(theta);
460 825 : theta += 2.0 * libMesh::pi / theta_res_square;
461 : }
462 : }
463 : std::array<Point, theta_res_triangle + 1> circle_points_triangle;
464 : {
465 : Real theta = 0;
466 858 : for (unsigned int i = 0; i < theta_res_triangle + 1; i++)
467 : {
468 825 : circle_points_triangle[i](0) = radius * std::cos(theta);
469 825 : circle_points_triangle[i](1) = radius * std::sin(theta);
470 825 : theta += 2.0 * libMesh::pi / theta_res_triangle;
471 : }
472 : }
473 : // Define "quadrant center" reference points. These will be the centers of
474 : // the 3 circles that represent the fuel pins. These centers are
475 : // offset a little bit so that in the final mesh, there is a tiny gap between
476 : // neighboring subchannel cells. That allows us to easily map a solution to
477 : // this detailed mesh with a nearest-neighbor search.
478 : const Real shrink_factor = 0.99999;
479 : // Quadrants are used only for the side and corner subchannels
480 33 : Real _duct_to_pin_gap =
481 33 : 0.5 * (_flat_to_flat - (_n_rings - 1) * _pitch * std::sqrt(3.0) - _pin_diameter);
482 : std::array<Point, 2> quadrant_centers_sides;
483 33 : quadrant_centers_sides[0] = Point(
484 33 : -_pitch * 0.5 * shrink_factor, -(_duct_to_pin_gap + _pin_diameter) * 0.5 * shrink_factor, 0);
485 33 : quadrant_centers_sides[1] = Point(
486 : _pitch * 0.5 * shrink_factor, -(_duct_to_pin_gap + _pin_diameter) * 0.5 * shrink_factor, 0);
487 : std::array<Point, 1> quadrant_centers_corner;
488 33 : quadrant_centers_corner[0] =
489 33 : Point(-(_duct_to_pin_gap + _pin_diameter) * 0.5 * std::sin(libMesh::pi / 6) * shrink_factor,
490 33 : -(_duct_to_pin_gap + _pin_diameter) * 0.5 * std::cos(libMesh::pi / 6) * shrink_factor,
491 : 0);
492 : // Triangles are used for all center subchannels
493 : std::array<Point, 3> triangle_centers;
494 33 : triangle_centers[0] = Point(0, _pitch * std::cos(libMesh::pi / 6) * 2 / 3 * shrink_factor, 0);
495 33 : triangle_centers[1] = Point(-_pitch * 0.5 * shrink_factor,
496 33 : -_pitch * std::cos(libMesh::pi / 6) * 1 / 3 * shrink_factor,
497 : 0);
498 33 : triangle_centers[2] = Point(
499 : _pitch * 0.5 * shrink_factor, -_pitch * std::cos(libMesh::pi / 6) * 1 / 3 * shrink_factor, 0);
500 :
501 : const unsigned int m_sixth = theta_res_triangle / 6;
502 : const unsigned int m_quarter = theta_res_square / 4;
503 : // Build an array of points that represent a cross section of a center subchannel
504 : // cell. The points are ordered in this fashion:
505 : // 3 1
506 : // 2
507 : // 0
508 : // 4 5 8 9
509 : // 6 7
510 561 : std::array<Point, points_per_center> center_points{};
511 33 : center_points[0] = Point(0, 0, 0);
512 : {
513 : unsigned int start;
514 132 : for (unsigned int i = 0; i < 3; i++)
515 : {
516 99 : if (i == 0)
517 : start = 5 * (m_sixth);
518 66 : if (i == 1)
519 : start = 1 * (m_sixth);
520 66 : if (i == 2)
521 : start = 3 * (m_sixth);
522 594 : for (unsigned int ii = 0; ii < points_per_sixth; ii++)
523 : {
524 495 : auto c_pt = circle_points_triangle[start - ii];
525 495 : center_points[i * points_per_sixth + ii + 1] = triangle_centers[i] + c_pt;
526 : }
527 : }
528 : }
529 :
530 : // Build an array of points that represent a cross section of a top left corner subchannel
531 : // cell. The points are ordered in this fashion (x direction towards point 4, y direction towards
532 : // point 6):
533 : // 5 6
534 : //
535 : // 0
536 : // 4 2 3
537 : // 1
538 330 : std::array<Point, points_per_corner> corner_points{};
539 33 : corner_points[0] = Point(0, 0, 0);
540 : {
541 198 : for (unsigned int ii = 0; ii < points_per_sixth; ii++)
542 : {
543 165 : auto c_pt = circle_points_triangle[1 * m_quarter - ii];
544 165 : corner_points[ii + 1] = quadrant_centers_corner[0] + c_pt;
545 : }
546 : Real side_short = (_duct_to_pin_gap + _pin_diameter) * 0.5;
547 33 : Real side_long = (2.0 * _duct_to_pin_gap + _pin_diameter) * 0.5;
548 33 : Real side_length = std::sqrt(std::pow(side_short, 2) + std::pow(side_long, 2) -
549 33 : 2 * side_short * side_long * std::cos(libMesh::pi / 6));
550 : Real angle =
551 : libMesh::pi - libMesh::pi / 3 -
552 33 : std::acos((-std::pow(side_long, 2) + std::pow(side_short, 2) + std::pow(side_length, 2)) /
553 33 : (2 * side_short * side_length));
554 33 : corner_points[points_per_sixth + 1] = Point(side_length * std::cos(angle) * shrink_factor,
555 33 : -side_length * std::sin(angle) * shrink_factor,
556 : 0);
557 33 : corner_points[points_per_sixth + 2] =
558 33 : Point(side_length * std::sin(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor *
559 : std::tan(libMesh::pi / 6),
560 33 : side_length * std::sin(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor,
561 : 0);
562 33 : corner_points[points_per_sixth + 3] =
563 33 : Point(-side_length * std::cos(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor,
564 : side_length * std::sin(libMesh::pi / 2 - angle - libMesh::pi / 6) * shrink_factor,
565 : 0);
566 : }
567 :
568 : // Build an array of points that represent a cross-section of a top side subchannel
569 : // cell. The points are ordered in this fashion:
570 : // 8 7
571 : //
572 : // 0
573 : // 1 2 5 6
574 : // 3 4
575 594 : std::array<Point, points_per_side> side_points{};
576 33 : side_points[0] = Point(0, 0, 0);
577 : {
578 264 : for (unsigned int ii = 0; ii < points_per_quadrant; ii++)
579 : {
580 231 : auto c_pt = circle_points_square[m_quarter - ii];
581 : side_points[ii + 1] = quadrant_centers_sides[0] + c_pt;
582 : }
583 264 : for (unsigned int ii = 0; ii < points_per_quadrant; ii++)
584 : {
585 231 : auto c_pt = circle_points_square[2 * m_quarter - ii];
586 231 : side_points[points_per_quadrant + ii + 1] = quadrant_centers_sides[1] + c_pt;
587 : }
588 33 : side_points[2 * points_per_quadrant + 1] =
589 33 : Point(_pitch * 0.5 * shrink_factor, 0.5 * _duct_to_pin_gap * shrink_factor, 0);
590 33 : side_points[2 * points_per_quadrant + 2] =
591 : Point(-_pitch * 0.5 * shrink_factor, 0.5 * _duct_to_pin_gap * shrink_factor, 0);
592 : }
593 :
594 33 : int point_counter = 0;
595 : unsigned int node_id = 0;
596 1635 : for (unsigned int i = 0; i < _n_channels; i++)
597 : {
598 : // Real offset_x = _flat_to_flat / 2.0;
599 : // Real offset_y = _flat_to_flat / 2.0;
600 :
601 1602 : if (getSubchannelType(i) == EChannelType::CENTER)
602 : {
603 34344 : for (auto z : _z_grid)
604 : {
605 : // Get height
606 : Point z0{0, 0, z};
607 :
608 : // Get suchannel position and assign to point
609 33372 : auto loc_position = getSubchannelPosition(i);
610 33372 : Point p0{loc_position[0], loc_position[1], 0};
611 :
612 : // Determine orientation of current subchannel
613 33372 : auto subchannel_pins = getSubChannelPins(i);
614 : Point subchannel_side =
615 33372 : getPinPosition(subchannel_pins[0]) + getPinPosition(subchannel_pins[1]);
616 : Point base_center_orientation = {0, -1};
617 :
618 : // Get rotation angle for current subchannel
619 : Real dot_prod = 0;
620 100116 : for (unsigned int lp = 0; lp < 2; lp++)
621 66744 : dot_prod += base_center_orientation(lp) * subchannel_side(lp);
622 : auto theta =
623 33372 : std::acos(dot_prod / (base_center_orientation.norm() * subchannel_side.norm()));
624 33372 : if (subchannel_side(0) < 0)
625 16686 : theta = 2.0 * libMesh::pi - theta;
626 :
627 : // Real distance_side = subchannel_side.norm();
628 : // Real distance_top = getPinPosition(subchannel_pins[2]).norm();
629 : // if (distance_top > distance_side)
630 : // theta += libMesh::pi * 0.0;
631 :
632 33372 : theta += _orientation_map[i];
633 :
634 33372 : theta = trunc((theta + (libMesh::pi / 6.0)) / (libMesh::pi / 3.0)) * libMesh::pi / 3.0;
635 :
636 33372 : if (_verbose)
637 : {
638 3528 : if (z == 0)
639 : {
640 168 : _console << "Subchannel Position: " << p0 << std::endl;
641 168 : auto pins = getSubChannelPins(i);
642 672 : for (auto r : pins)
643 504 : _console << r << " ";
644 168 : _console << std::endl;
645 168 : _console << "Theta: " << theta / libMesh::pi * 180. << std::endl;
646 168 : }
647 : }
648 :
649 : // Assigning points for center channels
650 567324 : for (unsigned int i = 0; i < points_per_center; i++)
651 : {
652 533952 : auto new_point = rotatePoint(center_points[i], theta) + p0 + z0;
653 533952 : if (_verbose)
654 : {
655 56448 : if (z == 0)
656 2688 : _console << i << " - " << new_point << std::endl;
657 : }
658 533952 : mesh_base->add_point(new_point, node_id++);
659 533952 : point_counter += 1;
660 : }
661 33372 : }
662 : }
663 630 : else if (getSubchannelType(i) == EChannelType::EDGE)
664 : {
665 15984 : for (auto z : _z_grid)
666 : {
667 : // Get height
668 : Point z0{0, 0, z};
669 :
670 : // Get suchannel position and assign to point
671 15552 : auto loc_position = getSubchannelPosition(i);
672 15552 : Point p0{loc_position[0], loc_position[1], 0};
673 :
674 : // Determine orientation of current subchannel
675 15552 : auto subchannel_pins = getSubChannelPins(i);
676 : Point subchannel_side =
677 15552 : getPinPosition(subchannel_pins[0]) + getPinPosition(subchannel_pins[1]);
678 : Point base_center_orientation = {0, 1};
679 :
680 : // Get rotation angle for current subchannel
681 : Real dot_prod = 0;
682 46656 : for (unsigned int lp = 0; lp < 2; lp++)
683 31104 : dot_prod += base_center_orientation(lp) * subchannel_side(lp);
684 : auto theta =
685 15552 : std::acos(dot_prod / (base_center_orientation.norm() * subchannel_side.norm()));
686 15552 : if (subchannel_side(0) > 0)
687 7776 : theta = 2. * libMesh::pi - theta;
688 15552 : theta = trunc((theta + (libMesh::pi / 6.0)) / (libMesh::pi / 3.0)) * libMesh::pi / 3.0;
689 :
690 15552 : if (_verbose)
691 : {
692 1764 : if (z == 0)
693 : {
694 84 : _console << "Subchannel Position: " << p0 << std::endl;
695 84 : auto pins = getSubChannelPins(i);
696 252 : for (auto r : pins)
697 168 : _console << r << " ";
698 84 : _console << std::endl;
699 84 : _console << "Theta: " << theta * 180 / libMesh::pi << std::endl;
700 84 : }
701 : }
702 :
703 : // Assigning points for center channels
704 279936 : for (unsigned int i = 0; i < points_per_side; i++)
705 : {
706 264384 : auto new_point = rotatePoint(side_points[i], theta) + p0 + z0;
707 264384 : if (_verbose)
708 : {
709 29988 : if (z == 0)
710 1428 : _console << i << " - " << new_point << std::endl;
711 : }
712 264384 : mesh_base->add_point(new_point, node_id++);
713 264384 : point_counter += 1;
714 : }
715 15552 : }
716 : }
717 : else // getSubchannelType(i) == EChannelType::CORNER
718 : {
719 7596 : for (auto z : _z_grid)
720 : {
721 : // Get height
722 : Point z0{0, 0, z};
723 :
724 : // Get suchannel position and assign to point
725 7398 : auto loc_position = getSubchannelPosition(i);
726 7398 : Point p0{loc_position[0], loc_position[1], 0};
727 :
728 : // Determine orientation of current subchannel
729 7398 : auto subchannel_pins = getSubChannelPins(i);
730 7398 : Point subchannel_side = getPinPosition(subchannel_pins[0]);
731 : Point base_center_orientation = {1, 1};
732 :
733 : // Get rotation angle for current subchannel
734 : Real dot_prod = 0;
735 22194 : for (unsigned int lp = 0; lp < 2; lp++)
736 14796 : dot_prod += base_center_orientation(lp) * subchannel_side(lp);
737 : auto theta =
738 7398 : std::acos(dot_prod / (base_center_orientation.norm() * subchannel_side.norm()));
739 7398 : if (subchannel_side(0) > 0)
740 3699 : theta = 2. * libMesh::pi - theta;
741 7398 : theta = trunc((theta + (libMesh::pi / 6.0)) / (libMesh::pi / 3.0)) * libMesh::pi / 3.0;
742 :
743 7398 : if (_verbose)
744 : {
745 882 : if (z == 0)
746 : {
747 42 : _console << "Subchannel Position: " << p0 << std::endl;
748 42 : auto pins = getSubChannelPins(i);
749 84 : for (auto r : pins)
750 42 : _console << r << " ";
751 42 : _console << std::endl;
752 42 : _console << "Theta: " << theta * 180 / libMesh::pi << std::endl;
753 42 : }
754 : }
755 :
756 : // Assigning points for center channels
757 73980 : for (unsigned int i = 0; i < points_per_corner; i++)
758 : {
759 66582 : auto new_point = rotatePoint(corner_points[i], theta) + p0 + z0;
760 66582 : if (_verbose)
761 : {
762 7938 : if (z == 0)
763 378 : _console << i << " - " << new_point << std::endl;
764 : }
765 66582 : mesh_base->add_point(new_point, node_id++);
766 66582 : point_counter += 1;
767 : }
768 7398 : }
769 : }
770 : } // i
771 33 : if (_verbose)
772 7 : _console << "Point counter: " << point_counter << std::endl;
773 :
774 33 : int element_counter = 0;
775 : unsigned int elem_id = 0;
776 : unsigned int number_of_corner = 0;
777 : unsigned int number_of_side = 0;
778 : unsigned int number_of_center = 0;
779 : unsigned int elems_per_channel = 0;
780 : unsigned int points_per_channel = 0;
781 1635 : for (unsigned int i = 0; i < _n_channels; i++)
782 : {
783 : auto subch_type = getSubchannelType(i);
784 1602 : if (subch_type == EChannelType::CORNER)
785 : {
786 198 : number_of_corner++;
787 : elems_per_channel = elems_per_corner;
788 : points_per_channel = points_per_corner;
789 198 : if (_verbose)
790 42 : _console << "Corner" << std::endl;
791 : }
792 1404 : else if (subch_type == EChannelType::EDGE)
793 : {
794 432 : number_of_side++;
795 : elems_per_channel = elems_per_side;
796 : points_per_channel = points_per_side;
797 432 : if (_verbose)
798 84 : _console << "Edge" << std::endl;
799 : }
800 972 : else if (subch_type == EChannelType::CENTER)
801 : {
802 972 : number_of_center++;
803 : elems_per_channel = elems_per_center;
804 : points_per_channel = points_per_center;
805 972 : if (_verbose)
806 168 : _console << "Center" << std::endl;
807 : }
808 56322 : for (unsigned int iz = 0; iz < _n_cells; iz++)
809 : {
810 54720 : unsigned int elapsed_points = number_of_corner * points_per_corner * (_n_cells + 1) +
811 54720 : number_of_side * points_per_side * (_n_cells + 1) +
812 54720 : number_of_center * points_per_center * (_n_cells + 1) -
813 54720 : points_per_channel * (_n_cells + 1);
814 : // index of the central node at base of cell
815 54720 : unsigned int indx1 = iz * points_per_channel + elapsed_points;
816 : // index of the central node at top of cell
817 54720 : unsigned int indx2 = (iz + 1) * points_per_channel + elapsed_points;
818 :
819 840240 : for (unsigned int i = 0; i < elems_per_channel; i++)
820 : {
821 785520 : Elem * elem = mesh_base->add_elem(std::make_unique<Prism6>());
822 785520 : elem->subdomain_id() = _subchannel_block_id;
823 785520 : elem->set_id(elem_id++);
824 :
825 785520 : if (_verbose)
826 84000 : _console << "Node 0: " << *mesh_base->node_ptr(indx1) << std::endl;
827 :
828 785520 : elem->set_node(0, mesh_base->node_ptr(indx1));
829 785520 : elem->set_node(1, mesh_base->node_ptr(indx1 + i + 1));
830 785520 : if (i != elems_per_channel - 1)
831 730800 : elem->set_node(2, mesh_base->node_ptr(indx1 + i + 2));
832 : else
833 54720 : elem->set_node(2, mesh_base->node_ptr(indx1 + 1));
834 :
835 785520 : elem->set_node(3, mesh_base->node_ptr(indx2));
836 785520 : elem->set_node(4, mesh_base->node_ptr(indx2 + i + 1));
837 785520 : if (i != elems_per_channel - 1)
838 730800 : elem->set_node(5, mesh_base->node_ptr(indx2 + i + 2));
839 : else
840 54720 : elem->set_node(5, mesh_base->node_ptr(indx2 + 1));
841 :
842 785520 : if (iz == 0)
843 23076 : boundary_info.add_side(elem, 0, 0);
844 785520 : if (iz == _n_cells - 1)
845 23076 : boundary_info.add_side(elem, 4, 1);
846 :
847 785520 : element_counter += 1;
848 : }
849 : }
850 : }
851 33 : if (_verbose)
852 7 : _console << "Element counter: " << element_counter << std::endl;
853 33 : boundary_info.sideset_name(0) = "inlet";
854 33 : boundary_info.sideset_name(1) = "outlet";
855 33 : mesh_base->subdomain_name(_subchannel_block_id) = "subchannel";
856 33 : if (n_pins > 0)
857 : {
858 33 : _elem_id = mesh_base->n_elem();
859 768 : for (auto & ctr : _pin_position)
860 735 : generatePin(mesh_base, ctr);
861 : }
862 : if (n_pins > 0)
863 33 : mesh_base->subdomain_name(_pin_block_id) = "fuel_pins";
864 33 : if (_verbose)
865 7 : _console << "Mesh assembly done" << std::endl;
866 33 : mesh_base->prepare_for_use();
867 :
868 33 : return mesh_base;
869 0 : }
870 :
871 : Point
872 864918 : SCMDetailedTriAssemblyMeshGenerator::rotatePoint(Point b, Real theta)
873 : {
874 864918 : std::vector<std::vector<Real>> A(3, std::vector<Real>(3));
875 :
876 864918 : A[0] = {std::cos(theta), -std::sin(theta), 0.0};
877 864918 : A[1] = {std::sin(theta), std::cos(theta), 0.0};
878 864918 : A[2] = {0.0, 0.0, 1.0};
879 :
880 : Point rotated_vector = Point(0.0, 0.0, 0.0);
881 3459672 : for (unsigned int i = 0; i < 3; i++)
882 10379016 : for (unsigned int j = 0; j < 3; j++)
883 7784262 : rotated_vector(i) += A[i][j] * b(j);
884 :
885 864918 : return rotated_vector;
886 864918 : }
|
