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 "SCMTriAssemblyMeshGenerator.h"
11 : #include "TriSubChannelMesh.h"
12 : #include <cmath>
13 : #include <memory>
14 : #include "libmesh/edge_edge2.h"
15 : #include "libmesh/unstructured_mesh.h"
16 :
17 : registerMooseObject("SubChannelApp", SCMTriAssemblyMeshGenerator);
18 : registerMooseObjectRenamed("SubChannelApp",
19 : SCMTriSubChannelMeshGenerator,
20 : "06/30/2027 24:00",
21 : SCMTriAssemblyMeshGenerator);
22 : registerMooseObjectRenamed("SubChannelApp",
23 : TriSubChannelMeshGenerator,
24 : "06/30/2027 24:00",
25 : SCMTriAssemblyMeshGenerator);
26 : registerMooseObjectRenamed("SubChannelApp",
27 : SCMTriPinMeshGenerator,
28 : "06/30/2027 24:00",
29 : SCMTriAssemblyMeshGenerator);
30 : registerMooseObjectRenamed("SubChannelApp",
31 : TriPinMeshGenerator,
32 : "06/30/2027 24:00",
33 : SCMTriAssemblyMeshGenerator);
34 :
35 : InputParameters
36 396 : SCMTriAssemblyMeshGenerator::validParams()
37 : {
38 396 : InputParameters params = MeshGenerator::validParams();
39 396 : params.addClassDescription(
40 : "Creates a mesh of 1D subchannels and 1D pins in a triangular lattice arrangement");
41 792 : params.addRequiredParam<unsigned int>("n_cells", "The number of cells in the axial direction");
42 792 : params.addRequiredParam<Real>("pitch", "Pitch [m]");
43 792 : params.addRequiredParam<Real>("pin_diameter", "Rod diameter [m]");
44 792 : params.addParam<Real>("unheated_length_entry", 0.0, "Unheated length at entry [m]");
45 792 : params.addRequiredParam<Real>("heated_length", "Heated length [m]");
46 792 : params.addParam<Real>("unheated_length_exit", 0.0, "Unheated length at exit [m]");
47 792 : params.addRequiredParam<unsigned int>(
48 : "nrings",
49 : "Number of fuel-pin rings per assembly, counting the center pin as the first ring [-]");
50 792 : params.addRequiredParam<Real>("flat_to_flat",
51 : "Flat to flat distance for the hexagonal assembly [m]");
52 792 : params.addRequiredParam<Real>("dwire", "Wire diameter [m]");
53 792 : params.addRequiredParam<Real>("hwire", "Wire lead length [m]");
54 792 : params.addParam<std::vector<Real>>(
55 : "spacer_z", {}, "Axial location of spacers/vanes/mixing vanes [m]");
56 792 : params.addParam<std::vector<Real>>(
57 : "spacer_k", {}, "K-loss coefficient of spacers/vanes/mixing vanes [-]");
58 792 : params.addParam<Real>("Kij", 0.5, "Lateral form loss coefficient [-]");
59 792 : params.addParam<std::vector<Real>>("z_blockage",
60 792 : std::vector<Real>({0.0, 0.0}),
61 : "axial location of blockage (inlet, outlet) [m]");
62 792 : params.addParam<std::vector<unsigned int>>("index_blockage",
63 792 : std::vector<unsigned int>({0}),
64 : "index of subchannels affected by blockage");
65 792 : params.addParam<std::vector<Real>>(
66 : "reduction_blockage",
67 792 : std::vector<Real>({1.0}),
68 : "Area reduction of subchannels affected by blockage (number to muliply the area)");
69 792 : params.addParam<std::vector<Real>>("k_blockage",
70 792 : std::vector<Real>({0.0}),
71 : "Form loss coefficient of subchannels affected by blockage");
72 792 : params.addParam<unsigned int>("block_id", 0, "Subchannel block id");
73 792 : params.deprecateParam("block_id", "subchannel_block_id", "07/01/2027");
74 792 : params.addParam<unsigned int>("pin_block_id", 1, "Fuel Pin block id");
75 396 : return params;
76 0 : }
77 :
78 198 : SCMTriAssemblyMeshGenerator::SCMTriAssemblyMeshGenerator(const InputParameters & params)
79 : : MeshGenerator(params),
80 198 : _unheated_length_entry(getParam<Real>("unheated_length_entry")),
81 396 : _heated_length(getParam<Real>("heated_length")),
82 396 : _unheated_length_exit(getParam<Real>("unheated_length_exit")),
83 396 : _subchannel_block_id(getParam<unsigned int>("subchannel_block_id")),
84 396 : _pin_block_id(getParam<unsigned int>("pin_block_id")),
85 396 : _spacer_z(getParam<std::vector<Real>>("spacer_z")),
86 396 : _spacer_k(getParam<std::vector<Real>>("spacer_k")),
87 396 : _z_blockage(getParam<std::vector<Real>>("z_blockage")),
88 396 : _index_blockage(getParam<std::vector<unsigned int>>("index_blockage")),
89 396 : _reduction_blockage(getParam<std::vector<Real>>("reduction_blockage")),
90 396 : _k_blockage(getParam<std::vector<Real>>("k_blockage")),
91 396 : _pitch(getParam<Real>("pitch")),
92 396 : _kij(getParam<Real>("Kij")),
93 396 : _pin_diameter(getParam<Real>("pin_diameter")),
94 396 : _n_cells(getParam<unsigned int>("n_cells")),
95 396 : _n_rings(getParam<unsigned int>("nrings")),
96 198 : _n_channels(0),
97 396 : _flat_to_flat(getParam<Real>("flat_to_flat")),
98 396 : _dwire(getParam<Real>("dwire")),
99 396 : _hwire(getParam<Real>("hwire")),
100 198 : _duct_to_pin_gap(0.5 *
101 198 : (_flat_to_flat - (_n_rings - 1) * _pitch * std::sqrt(3.0) - _pin_diameter)),
102 198 : _npins(0),
103 198 : _n_gaps(0)
104 : {
105 198 : const Real total_length = _unheated_length_entry + _heated_length + _unheated_length_exit;
106 :
107 198 : if (_n_rings < 2)
108 0 : paramError("nrings",
109 : "'nrings' must be at least 2. In this mesh generator, the center pin counts as "
110 : "the first ring, so a 7-pin bundle uses nrings = 2.");
111 :
112 198 : if (_n_cells == 0)
113 0 : paramError("n_cells", "The number of axial cells must be greater than zero");
114 :
115 198 : if (total_length <= 0.0)
116 0 : mooseError("Total bundle length must be greater than zero");
117 :
118 198 : if (_spacer_z.size() != _spacer_k.size())
119 0 : mooseError("Size of vector spacer_z should be equal to size of vector spacer_k");
120 :
121 382 : for (const auto spacer_z : _spacer_z)
122 184 : if (spacer_z < 0.0 || spacer_z > total_length)
123 0 : paramError("spacer_z", "Location of spacers should be between zero and total bundle length");
124 :
125 198 : if (_z_blockage.size() != 2)
126 0 : paramError("z_blockage", "Size of vector z_blockage must be 2");
127 :
128 198 : if (_z_blockage.front() > _z_blockage.back())
129 0 : paramError("z_blockage", "z_blockage inlet location must not exceed outlet location");
130 :
131 198 : if (*max_element(_reduction_blockage.begin(), _reduction_blockage.end()) > 1)
132 0 : paramError("reduction_blockage",
133 : "The area reduction of the blocked subchannels cannot be more than 1");
134 :
135 198 : if ((_index_blockage.size() != _reduction_blockage.size()) ||
136 396 : (_index_blockage.size() != _k_blockage.size()) ||
137 : (_reduction_blockage.size() != _k_blockage.size()))
138 0 : mooseError("Size of vectors: index_blockage, reduction_blockage, k_blockage, must be equal "
139 : "to eachother");
140 :
141 198 : SubChannelMesh::generateZGrid(
142 198 : _unheated_length_entry, _heated_length, _unheated_length_exit, _n_cells, _z_grid);
143 :
144 : // Defining the total length from 3 axial sections
145 : Real L = total_length;
146 :
147 : // Defining the position of the spacer grid in the numerical solution array
148 : std::vector<int> spacer_cell;
149 382 : for (const auto & elem : _spacer_z)
150 184 : spacer_cell.emplace_back(std::round(elem * _n_cells / L));
151 :
152 : // Defining the array for axial resistances
153 : std::vector<Real> kgrid;
154 198 : kgrid.resize(_n_cells + 1, 0.0);
155 :
156 : // Summing the spacer resistance to the grid resistance array
157 382 : for (unsigned int index = 0; index < spacer_cell.size(); index++)
158 184 : kgrid[spacer_cell[index]] += _spacer_k[index];
159 :
160 : // compute the hex mesh variables
161 : // -------------------------------------------
162 : // x coordinate for the first position
163 : Real x0 = 0.0;
164 : // y coordinate for the first position
165 : Real y0 = 0.0;
166 : // x coordinate for the second position
167 : Real x1 = 0.0;
168 : // y coordinate for the second position dummy variable
169 : Real y1 = 0.0;
170 : // dummy variable
171 : Real a1 = 0.0;
172 : // dummy variable
173 : Real a2 = 0.0;
174 : // average x coordinate
175 : Real avg_coor_x = 0.0;
176 : // average y coordinate
177 : Real avg_coor_y = 0.0;
178 : // distance between two points
179 : Real dist = 0.0;
180 : // distance between two points
181 : Real dist0 = 0.0;
182 : // integer counter
183 198 : unsigned int kgap = 0;
184 : // dummy integer
185 : unsigned int icorner = 0;
186 : // used to defined global direction of the cross_flow_map coefficients for each subchannel and gap
187 198 : const Real positive_flow = 1.0;
188 : // used to defined global direction of the cross_flow_map coefficients for each subchannel and gap
189 : const Real negative_flow = -1.0;
190 : // the indicator used while setting _gap_to_chan_map array
191 : std::vector<std::pair<unsigned int, unsigned int>> gap_fill;
192 198 : TriSubChannelMesh::pinPositions(_pin_position, _n_rings, _pitch, Point(0, 0));
193 198 : _npins = _pin_position.size();
194 : // assign the pins to the corresponding rings
195 : unsigned int k = 0; // initialize the fuel Pin counter index
196 198 : _pins_in_rings.resize(_n_rings);
197 198 : _pins_in_rings[0].push_back(k++);
198 755 : for (unsigned int i = 1; i < _n_rings; i++)
199 8195 : for (unsigned int j = 0; j < i * 6; j++)
200 7638 : _pins_in_rings[i].push_back(k++);
201 : // Given the number of pins and number of fuel Pin rings, the number of subchannels can be
202 : // computed as follows:
203 : unsigned int chancount = 0.0;
204 : // Summing internal channels
205 755 : for (unsigned int j = 0; j < _n_rings - 1; j++)
206 557 : chancount += j * 6;
207 : // Adding external channels to the total count
208 198 : _n_channels = chancount + _npins - 1 + (_n_rings - 1) * 6 + 6;
209 :
210 198 : if (*max_element(_index_blockage.begin(), _index_blockage.end()) > (_n_channels - 1))
211 0 : paramError("index_blockage",
212 : "The index of the blocked subchannel cannot be more than the max index of the "
213 : "subchannels");
214 :
215 198 : if ((_index_blockage.size() > _n_channels) || (_reduction_blockage.size() > _n_channels) ||
216 : (_k_blockage.size() > _n_channels))
217 0 : mooseError("Size of vectors: index_blockage, reduction_blockage, k_blockage, cannot be more "
218 : "than the total number of subchannels");
219 :
220 : // Defining the 2D array for axial resistances
221 198 : _k_grid.resize(_n_channels, std::vector<Real>(_n_cells + 1));
222 16662 : for (unsigned int i = 0; i < _n_channels; i++)
223 16464 : _k_grid[i] = kgrid;
224 :
225 : // Add blockage resistance to the 2D grid resistane array
226 198 : Real dz = L / _n_cells;
227 4887 : for (unsigned int i = 0; i < _n_cells + 1; i++)
228 : {
229 4689 : if ((dz * i >= _z_blockage.front() && dz * i <= _z_blockage.back()))
230 : {
231 : unsigned int index(0);
232 613 : for (const auto & i_ch : _index_blockage)
233 : {
234 415 : _k_grid[i_ch][i] += _k_blockage[index];
235 415 : index++;
236 : }
237 : }
238 : }
239 :
240 198 : _chan_to_pin_map.resize(_n_channels);
241 198 : _pin_to_chan_map.resize(_npins);
242 198 : _subch_type.resize(_n_channels);
243 198 : _n_gaps = _n_channels + _npins - 1; /// initial assignment
244 198 : _gap_to_chan_map.resize(_n_gaps);
245 198 : _gap_to_pin_map.resize(_n_gaps);
246 198 : gap_fill.resize(_n_gaps);
247 198 : _chan_to_gap_map.resize(_n_channels);
248 198 : _gap_pairs_sf.resize(_n_channels);
249 198 : _chan_pairs_sf.resize(_n_channels);
250 198 : _gij_map.resize(_n_cells + 1);
251 198 : _sign_id_crossflow_map.resize(_n_channels);
252 198 : _gap_type.resize(_n_gaps);
253 198 : _subchannel_position.resize(_n_channels);
254 :
255 16662 : for (unsigned int i = 0; i < _n_channels; i++)
256 : {
257 16464 : _chan_to_pin_map[i].reserve(3);
258 16464 : _chan_to_gap_map[i].reserve(3);
259 16464 : _sign_id_crossflow_map[i].reserve(3);
260 16464 : _subchannel_position[i].reserve(3);
261 65856 : for (unsigned int j = 0; j < 3; j++)
262 : {
263 49392 : _sign_id_crossflow_map.at(i).push_back(positive_flow);
264 49392 : _subchannel_position.at(i).push_back(0.0);
265 : }
266 : }
267 :
268 4887 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
269 : {
270 4689 : _gij_map[iz].reserve(_n_gaps);
271 : }
272 :
273 8034 : for (unsigned int i = 0; i < _npins; i++)
274 7836 : _pin_to_chan_map[i].reserve(6);
275 :
276 : // create the subchannels
277 : k = 0; // initialize the subchannel counter index
278 : kgap = 0;
279 : // for each ring we trace the subchannels by pairing up to neighbor pins and looking for the third
280 : // Pin at inner or outer ring compared to the current ring.
281 755 : for (unsigned int i = 1; i < _n_rings; i++)
282 : {
283 : // find the closest Pin at back ring
284 8195 : for (unsigned int j = 0; j < _pins_in_rings[i].size(); j++)
285 : {
286 7638 : if (j == _pins_in_rings[i].size() - 1)
287 : {
288 557 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
289 557 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][0]);
290 557 : avg_coor_x =
291 557 : 0.5 * (_pin_position[_pins_in_rings[i][j]](0) + _pin_position[_pins_in_rings[i][0]](0));
292 557 : avg_coor_y =
293 557 : 0.5 * (_pin_position[_pins_in_rings[i][j]](1) + _pin_position[_pins_in_rings[i][0]](1));
294 557 : _gap_to_pin_map[kgap].first = _pins_in_rings[i][0];
295 557 : _gap_to_pin_map[kgap].second = _pins_in_rings[i][j];
296 557 : _gap_type[kgap] = EChannelType::CENTER;
297 557 : kgap = kgap + 1;
298 : }
299 : else
300 : {
301 7081 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
302 7081 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j + 1]);
303 7081 : avg_coor_x = 0.5 * (_pin_position[_pins_in_rings[i][j]](0) +
304 7081 : _pin_position[_pins_in_rings[i][j + 1]](0));
305 7081 : avg_coor_y = 0.5 * (_pin_position[_pins_in_rings[i][j]](1) +
306 7081 : _pin_position[_pins_in_rings[i][j + 1]](1));
307 7081 : _gap_to_pin_map[kgap].first = _pins_in_rings[i][j];
308 7081 : _gap_to_pin_map[kgap].second = _pins_in_rings[i][j + 1];
309 7081 : _gap_type[kgap] = EChannelType::CENTER;
310 7081 : kgap = kgap + 1;
311 : }
312 :
313 : dist0 = 1.0e+5;
314 :
315 7638 : _chan_to_pin_map[k].push_back(_pins_in_rings[i - 1][0]);
316 : unsigned int l0 = 0;
317 :
318 109986 : for (unsigned int l = 0; l < _pins_in_rings[i - 1].size(); l++)
319 : {
320 102348 : dist = std::sqrt(pow(_pin_position[_pins_in_rings[i - 1][l]](0) - avg_coor_x, 2) +
321 102348 : pow(_pin_position[_pins_in_rings[i - 1][l]](1) - avg_coor_y, 2));
322 :
323 102348 : if (dist < dist0)
324 : {
325 34785 : _chan_to_pin_map[k][2] = _pins_in_rings[i - 1][l];
326 : l0 = l;
327 : dist0 = dist;
328 : } // if
329 : } // l
330 :
331 7638 : _gap_to_pin_map[kgap].first = _pins_in_rings[i][j];
332 7638 : _gap_to_pin_map[kgap].second = _pins_in_rings[i - 1][l0];
333 7638 : _gap_type[kgap] = EChannelType::CENTER;
334 7638 : kgap = kgap + 1;
335 7638 : _subch_type[k] = EChannelType::CENTER;
336 7638 : k = k + 1;
337 : } // for j
338 :
339 : // find the closest Pin at front ring
340 8195 : for (unsigned int j = 0; j < _pins_in_rings[i].size(); j++)
341 : {
342 7638 : if (j == _pins_in_rings[i].size() - 1)
343 : {
344 557 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
345 557 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][0]);
346 557 : avg_coor_x =
347 557 : 0.5 * (_pin_position[_pins_in_rings[i][j]](0) + _pin_position[_pins_in_rings[i][0]](0));
348 557 : avg_coor_y =
349 557 : 0.5 * (_pin_position[_pins_in_rings[i][j]](1) + _pin_position[_pins_in_rings[i][0]](1));
350 : }
351 : else
352 : {
353 7081 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
354 7081 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j + 1]);
355 7081 : avg_coor_x = 0.5 * (_pin_position[_pins_in_rings[i][j]](0) +
356 7081 : _pin_position[_pins_in_rings[i][j + 1]](0));
357 7081 : avg_coor_y = 0.5 * (_pin_position[_pins_in_rings[i][j]](1) +
358 7081 : _pin_position[_pins_in_rings[i][j + 1]](1));
359 : }
360 :
361 : // if the outermost ring, set the edge subchannels first... then the corner subchannels
362 7638 : if (i == _n_rings - 1)
363 : {
364 : // add edges
365 3342 : _subch_type[k] = EChannelType::EDGE; // an edge subchannel is created
366 3342 : _gap_to_pin_map[kgap].first = _pins_in_rings[i][j];
367 3342 : _gap_to_pin_map[kgap].second = _pins_in_rings[i][j];
368 3342 : _gap_type[kgap] = EChannelType::EDGE;
369 3342 : _chan_to_gap_map[k].push_back(kgap);
370 3342 : kgap = kgap + 1;
371 3342 : k = k + 1;
372 :
373 3342 : if (j % i == 0)
374 : {
375 : // generate a corner subchannel, generate the additional gap and fix chan_to_gap_map
376 1188 : _gap_to_pin_map[kgap].first = _pins_in_rings[i][j];
377 1188 : _gap_to_pin_map[kgap].second = _pins_in_rings[i][j];
378 1188 : _gap_type[kgap] = EChannelType::CORNER;
379 :
380 : // corner subchannel
381 1188 : _chan_to_pin_map[k].push_back(_pins_in_rings[i][j]);
382 1188 : _chan_to_gap_map[k].push_back(kgap - 1);
383 1188 : _chan_to_gap_map[k].push_back(kgap);
384 1188 : _subch_type[k] = EChannelType::CORNER;
385 :
386 1188 : kgap = kgap + 1;
387 1188 : k = k + 1;
388 : }
389 : // if not the outer most ring
390 : }
391 : else
392 : {
393 : dist0 = 1.0e+5;
394 : unsigned int l0 = 0;
395 4296 : _chan_to_pin_map[k].push_back(_pins_in_rings[i + 1][0]);
396 105456 : for (unsigned int l = 0; l < _pins_in_rings[i + 1].size(); l++)
397 : {
398 101160 : dist = std::sqrt(pow(_pin_position[_pins_in_rings[i + 1][l]](0) - avg_coor_x, 2) +
399 101160 : pow(_pin_position[_pins_in_rings[i + 1][l]](1) - avg_coor_y, 2));
400 101160 : if (dist < dist0)
401 : {
402 31988 : _chan_to_pin_map[k][2] = _pins_in_rings[i + 1][l];
403 : dist0 = dist;
404 : l0 = l;
405 : } // if
406 : } // l
407 :
408 4296 : _gap_to_pin_map[kgap].first = _pins_in_rings[i][j];
409 4296 : _gap_to_pin_map[kgap].second = _pins_in_rings[i + 1][l0];
410 4296 : _gap_type[kgap] = EChannelType::CENTER;
411 4296 : kgap = kgap + 1;
412 4296 : _subch_type[k] = EChannelType::CENTER;
413 4296 : k = k + 1;
414 : } // if
415 : } // for j
416 : } // for i
417 :
418 : // Constructing pins to channels mao
419 8034 : for (unsigned int loc_rod = 0; loc_rod < _npins; loc_rod++)
420 : {
421 1315656 : for (unsigned int i = 0; i < _n_channels; i++)
422 : {
423 : bool rod_in_sc = false;
424 4936338 : for (unsigned int j : _chan_to_pin_map[i])
425 : {
426 3628518 : if (j == loc_rod)
427 : rod_in_sc = true;
428 : }
429 1307820 : if (rod_in_sc)
430 : {
431 43674 : _pin_to_chan_map[loc_rod].push_back(i);
432 : }
433 : }
434 : }
435 :
436 : /**
437 : * Build channel-gap connectivity from the channel-pin and gap-pin maps.
438 : *
439 : * Center channels are bounded by three center gaps connecting their three pin pairs. Edge
440 : * channels are bounded by one center gap between their two pins and two perimeter gaps along the
441 : * duct. Corner channels are bounded by the two perimeter gaps that meet at the corner pin.
442 : *
443 : * For a two-ring assembly every outer-ring pin is a corner pin. Consequently an edge channel can
444 : * find a corner channel at both of its endpoint pins. The edge still owns only three local gap
445 : * entries, matching _sign_id_crossflow_map and the reverse _gap_to_chan_map construction below,
446 : * so the second corner lookup is skipped once those three entries are already present.
447 : */
448 16662 : for (unsigned int i = 0; i < _n_channels; i++)
449 : {
450 16464 : if (_subch_type[i] == EChannelType::CENTER)
451 : {
452 3227418 : for (unsigned int j = 0; j < _n_gaps; j++)
453 : {
454 3215484 : if (_gap_type[j] == EChannelType::CENTER)
455 : {
456 2852808 : if (((_chan_to_pin_map[i][0] == _gap_to_pin_map[j].first) &&
457 2820348 : (_chan_to_pin_map[i][1] == _gap_to_pin_map[j].second)) ||
458 2809330 : ((_chan_to_pin_map[i][0] == _gap_to_pin_map[j].second) &&
459 31272 : (_chan_to_pin_map[i][1] == _gap_to_pin_map[j].first)))
460 : {
461 11934 : _chan_to_gap_map[i].push_back(j);
462 : }
463 :
464 2852808 : if (((_chan_to_pin_map[i][0] == _gap_to_pin_map[j].first) &&
465 2820348 : (_chan_to_pin_map[i][2] == _gap_to_pin_map[j].second)) ||
466 2808414 : ((_chan_to_pin_map[i][0] == _gap_to_pin_map[j].second) &&
467 31272 : (_chan_to_pin_map[i][2] == _gap_to_pin_map[j].first)))
468 : {
469 11934 : _chan_to_gap_map[i].push_back(j);
470 : }
471 :
472 2852808 : if (((_chan_to_pin_map[i][1] == _gap_to_pin_map[j].first) &&
473 2820348 : (_chan_to_pin_map[i][2] == _gap_to_pin_map[j].second)) ||
474 2817006 : ((_chan_to_pin_map[i][1] == _gap_to_pin_map[j].second) &&
475 31272 : (_chan_to_pin_map[i][2] == _gap_to_pin_map[j].first)))
476 : {
477 11934 : _chan_to_gap_map[i].push_back(j);
478 : }
479 : }
480 : } // for j
481 : }
482 4530 : else if (_subch_type[i] == EChannelType::EDGE)
483 : {
484 616098 : for (unsigned int j = 0; j < _n_gaps; j++)
485 : {
486 612756 : if (_gap_type[j] == EChannelType::CENTER)
487 : {
488 527784 : if (((_chan_to_pin_map[i][0] == _gap_to_pin_map[j].first) &&
489 521100 : (_chan_to_pin_map[i][1] == _gap_to_pin_map[j].second)) ||
490 517956 : ((_chan_to_pin_map[i][0] == _gap_to_pin_map[j].second) &&
491 5496 : (_chan_to_pin_map[i][1] == _gap_to_pin_map[j].first)))
492 : {
493 3342 : _chan_to_gap_map[i].push_back(j);
494 : }
495 : }
496 : }
497 :
498 : icorner = 0;
499 403863 : for (unsigned int k = 0; k < _n_channels; k++)
500 : {
501 401709 : if (_subch_type[k] == EChannelType::CORNER &&
502 17082 : _chan_to_pin_map[i][1] == _chan_to_pin_map[k][0])
503 : {
504 1188 : _chan_to_gap_map[i].push_back(_chan_to_gap_map[k][1]);
505 : icorner = 1;
506 : break;
507 : } // if
508 : } // for
509 :
510 : // Check whether the edge channel's first pin is also a corner pin. This second corner lookup
511 : // is only needed while the edge channel still needs another perimeter gap.
512 3342 : if (_chan_to_gap_map[i].size() < 3)
513 : {
514 304395 : for (unsigned int k = 0; k < _n_channels; k++)
515 : {
516 303345 : if (_subch_type[k] == EChannelType::CORNER &&
517 10164 : _chan_to_pin_map[i][0] == _chan_to_pin_map[k][0])
518 : {
519 1104 : _chan_to_gap_map[i].push_back(_chan_to_gap_map[k][1] + 1);
520 : icorner = 1;
521 : break;
522 : }
523 : }
524 : }
525 :
526 2238 : if (icorner == 0)
527 : {
528 1050 : _chan_to_gap_map[i].push_back(_chan_to_gap_map[i][0] + 1);
529 : }
530 : }
531 : }
532 :
533 : // find gap_to_chan_map pair
534 24300 : for (unsigned int j = 0; j < _n_gaps; j++)
535 : {
536 3996954 : for (unsigned int i = 0; i < _n_channels; i++)
537 : {
538 3972852 : if (_subch_type[i] == EChannelType::CENTER || _subch_type[i] == EChannelType::EDGE)
539 : {
540 3828240 : if ((j == _chan_to_gap_map[i][0]) || (j == _chan_to_gap_map[i][1]) ||
541 3797688 : (j == _chan_to_gap_map[i][2]))
542 : {
543 45828 : if (_gap_to_chan_map[j].first == 0 && gap_fill[j].first == 0)
544 : {
545 23904 : _gap_to_chan_map[j].first = i;
546 23904 : gap_fill[j].first = 1;
547 : }
548 21924 : else if (_gap_to_chan_map[j].second == 0 && gap_fill[j].second == 0)
549 : {
550 21924 : _gap_to_chan_map[j].second = i;
551 21924 : gap_fill[j].second = 1;
552 : }
553 : else
554 : {
555 : }
556 : }
557 : }
558 144612 : else if (_subch_type[i] == EChannelType::CORNER)
559 : {
560 144612 : if ((j == _chan_to_gap_map[i][0]) || (j == _chan_to_gap_map[i][1]))
561 : {
562 2376 : if (_gap_to_chan_map[j].first == 0 && gap_fill[j].first == 0)
563 : {
564 198 : _gap_to_chan_map[j].first = i;
565 198 : gap_fill[j].first = 1;
566 : }
567 2178 : else if (_gap_to_chan_map[j].second == 0 && gap_fill[j].second == 0)
568 : {
569 2178 : _gap_to_chan_map[j].second = i;
570 2178 : gap_fill[j].second = 1;
571 : }
572 : else
573 : {
574 : }
575 : }
576 : }
577 : } // i
578 : } // j
579 :
580 16662 : for (unsigned int k = 0; k < _n_channels; k++)
581 : {
582 16464 : if (_subch_type[k] == EChannelType::EDGE)
583 : {
584 3342 : _gap_pairs_sf[k].first = _chan_to_gap_map[k][0];
585 3342 : _gap_pairs_sf[k].second = _chan_to_gap_map[k][2];
586 : auto k1 = _gap_pairs_sf[k].first;
587 : auto k2 = _gap_pairs_sf[k].second;
588 3342 : if (_gap_to_chan_map[k1].first == k)
589 : {
590 1188 : _chan_pairs_sf[k].first = _gap_to_chan_map[k1].second;
591 : }
592 : else
593 : {
594 2154 : _chan_pairs_sf[k].first = _gap_to_chan_map[k1].first;
595 : }
596 :
597 3342 : if (_gap_to_chan_map[k2].first == k)
598 : {
599 3144 : _chan_pairs_sf[k].second = _gap_to_chan_map[k2].second;
600 : }
601 : else
602 : {
603 198 : _chan_pairs_sf[k].second = _gap_to_chan_map[k2].first;
604 : }
605 : }
606 13122 : else if (_subch_type[k] == EChannelType::CORNER)
607 : {
608 1188 : _gap_pairs_sf[k].first = _chan_to_gap_map[k][1];
609 1188 : _gap_pairs_sf[k].second = _chan_to_gap_map[k][0];
610 :
611 : auto k1 = _gap_pairs_sf[k].first;
612 : auto k2 = _gap_pairs_sf[k].second;
613 :
614 1188 : if (_gap_to_chan_map[k1].first == k)
615 : {
616 198 : _chan_pairs_sf[k].first = _gap_to_chan_map[k1].second;
617 : }
618 : else
619 : {
620 990 : _chan_pairs_sf[k].first = _gap_to_chan_map[k1].first;
621 : }
622 :
623 1188 : if (_gap_to_chan_map[k2].first == k)
624 : {
625 0 : _chan_pairs_sf[k].second = _gap_to_chan_map[k2].second;
626 : }
627 : else
628 : {
629 1188 : _chan_pairs_sf[k].second = _gap_to_chan_map[k2].first;
630 : }
631 : }
632 : }
633 :
634 : // set the _gij_map
635 4887 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
636 : {
637 597267 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
638 : {
639 592578 : if (_gap_type[i_gap] == EChannelType::CENTER)
640 : {
641 483834 : _gij_map[iz].push_back(_pitch - _pin_diameter);
642 : }
643 108744 : else if (_gap_type[i_gap] == EChannelType::EDGE || _gap_type[i_gap] == EChannelType::CORNER)
644 : {
645 108744 : _gij_map[iz].push_back(_duct_to_pin_gap);
646 : }
647 : }
648 : }
649 :
650 16662 : for (unsigned int i = 0; i < _n_channels; i++)
651 : {
652 16464 : if (_subch_type[i] == EChannelType::CENTER || _subch_type[i] == EChannelType::EDGE)
653 : {
654 61104 : for (unsigned int k = 0; k < 3; k++)
655 : {
656 11530548 : for (unsigned int j = 0; j < _n_gaps; j++)
657 : {
658 11484720 : if (_chan_to_gap_map[i][k] == j && i == _gap_to_chan_map[j].first)
659 : {
660 23904 : if (i > _gap_to_chan_map[j].second)
661 : {
662 0 : _sign_id_crossflow_map[i][k] = negative_flow;
663 : }
664 : else
665 : {
666 23904 : _sign_id_crossflow_map[i][k] = positive_flow;
667 : }
668 : }
669 11460816 : else if (_chan_to_gap_map[i][k] == j && i == _gap_to_chan_map[j].second)
670 : {
671 21924 : if (i > _gap_to_chan_map[j].first)
672 : {
673 21924 : _sign_id_crossflow_map[i][k] = negative_flow;
674 : }
675 : else
676 : {
677 0 : _sign_id_crossflow_map[i][k] = positive_flow;
678 : }
679 : }
680 : } // j
681 : } // k
682 : }
683 1188 : else if (_subch_type[i] == EChannelType::CORNER)
684 : {
685 3564 : for (unsigned int k = 0; k < 2; k++)
686 : {
687 291600 : for (unsigned int j = 0; j < _n_gaps; j++)
688 : {
689 289224 : if (_chan_to_gap_map[i][k] == j && i == _gap_to_chan_map[j].first)
690 : {
691 198 : if (i > _gap_to_chan_map[j].second)
692 : {
693 0 : _sign_id_crossflow_map[i][k] = negative_flow;
694 : }
695 : else
696 : {
697 198 : _sign_id_crossflow_map[i][k] = positive_flow;
698 : }
699 : }
700 289026 : else if (_chan_to_gap_map[i][k] == j && i == _gap_to_chan_map[j].second)
701 : {
702 2178 : if (i > _gap_to_chan_map[j].first)
703 : {
704 2178 : _sign_id_crossflow_map[i][k] = negative_flow;
705 : }
706 : else
707 : {
708 0 : _sign_id_crossflow_map[i][k] = positive_flow;
709 : }
710 : }
711 : } // j
712 : } // k
713 : } // subch_type =2
714 : } // i
715 :
716 : // set the subchannel positions
717 16662 : for (unsigned int i = 0; i < _n_channels; i++)
718 : {
719 16464 : if (_subch_type[i] == EChannelType::CENTER)
720 : {
721 11934 : _subchannel_position[i][0] =
722 11934 : (_pin_position[_chan_to_pin_map[i][0]](0) + _pin_position[_chan_to_pin_map[i][1]](0) +
723 11934 : _pin_position[_chan_to_pin_map[i][2]](0)) /
724 : 3.0;
725 11934 : _subchannel_position[i][1] =
726 11934 : (_pin_position[_chan_to_pin_map[i][0]](1) + _pin_position[_chan_to_pin_map[i][1]](1) +
727 11934 : _pin_position[_chan_to_pin_map[i][2]](1)) /
728 : 3.0;
729 : }
730 4530 : else if (_subch_type[i] == EChannelType::EDGE)
731 : {
732 418530 : for (unsigned int j = 0; j < _n_channels; j++)
733 : {
734 415188 : if (_subch_type[j] == EChannelType::CENTER &&
735 323532 : ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][0] &&
736 3342 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][1]) ||
737 320190 : (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][1] &&
738 3342 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][0])))
739 : {
740 3342 : x0 = _pin_position[_chan_to_pin_map[j][2]](0);
741 3342 : y0 = _pin_position[_chan_to_pin_map[j][2]](1);
742 : }
743 411846 : else if (_subch_type[j] == EChannelType::CENTER &&
744 320190 : ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][0] &&
745 0 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][2]) ||
746 320190 : (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][2] &&
747 2154 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][0])))
748 : {
749 0 : x0 = _pin_position[_chan_to_pin_map[j][1]](0);
750 0 : y0 = _pin_position[_chan_to_pin_map[j][1]](1);
751 : }
752 411846 : else if (_subch_type[j] == EChannelType::CENTER &&
753 320190 : ((_chan_to_pin_map[i][0] == _chan_to_pin_map[j][1] &&
754 3342 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][2]) ||
755 320190 : (_chan_to_pin_map[i][0] == _chan_to_pin_map[j][2] &&
756 2154 : _chan_to_pin_map[i][1] == _chan_to_pin_map[j][1])))
757 : {
758 0 : x0 = _pin_position[_chan_to_pin_map[j][0]](0);
759 0 : y0 = _pin_position[_chan_to_pin_map[j][0]](1);
760 : }
761 415188 : x1 = 0.5 *
762 415188 : (_pin_position[_chan_to_pin_map[i][0]](0) + _pin_position[_chan_to_pin_map[i][1]](0));
763 415188 : y1 = 0.5 *
764 415188 : (_pin_position[_chan_to_pin_map[i][0]](1) + _pin_position[_chan_to_pin_map[i][1]](1));
765 415188 : a1 = _pin_diameter / 2.0 + _duct_to_pin_gap / 2.0;
766 415188 : a2 = std::sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) + a1;
767 415188 : _subchannel_position[i][0] = (a2 * x1 - a1 * x0) / (a2 - a1);
768 415188 : _subchannel_position[i][1] = (a2 * y1 - a1 * y0) / (a2 - a1);
769 : } // j
770 : }
771 1188 : else if (_subch_type[i] == EChannelType::CORNER)
772 : {
773 1188 : x0 = _pin_position[0](0);
774 1188 : y0 = _pin_position[0](1);
775 1188 : x1 = _pin_position[_chan_to_pin_map[i][0]](0);
776 1188 : y1 = _pin_position[_chan_to_pin_map[i][0]](1);
777 1188 : a1 = _pin_diameter / 2.0 + _duct_to_pin_gap / 2.0;
778 1188 : a2 = std::sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) + a1;
779 1188 : _subchannel_position[i][0] = (a2 * x1 - a1 * x0) / (a2 - a1);
780 1188 : _subchannel_position[i][1] = (a2 * y1 - a1 * y0) / (a2 - a1);
781 : }
782 : }
783 :
784 : // Reduce reserved memory in the channel-to-gap map.
785 16662 : for (auto & gap : _chan_to_gap_map)
786 : {
787 : gap.shrink_to_fit();
788 : }
789 198 : }
790 :
791 : void
792 198 : SCMTriAssemblyMeshGenerator::buildPinMesh(MeshBase & mesh_base)
793 : {
794 198 : if (_npins == 0)
795 : return;
796 :
797 198 : mesh_base.reserve_elem(_n_cells * _npins);
798 198 : mesh_base.reserve_nodes((_n_cells + 1) * _npins);
799 :
800 198 : _pin_nodes.clear();
801 198 : _pin_nodes.resize(_npins);
802 :
803 : // Defining the extent of the subchannel mesh to append the pin mesh to the current subchannel
804 : // mesh.
805 198 : const unsigned int node_sub = mesh_base.n_nodes();
806 198 : const unsigned int elem_sub = mesh_base.n_elem();
807 :
808 : // Add the points in the shape of a rectilinear grid. The grid is regular on the xy-plane at the
809 : // triangular lattice pin positions. The grid along z is also regular. Store pointers in the
810 : // _pin_nodes array so we can keep track of which points are in which pins.
811 : unsigned int node_id = node_sub;
812 8034 : for (unsigned int i = 0; i < _npins; i++)
813 : {
814 7836 : _pin_nodes[i].reserve(_n_cells + 1);
815 200673 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
816 192837 : _pin_nodes[i].push_back(mesh_base.add_point(
817 385674 : Point(_pin_position[i](0), _pin_position[i](1), _z_grid[iz]), node_id++));
818 : }
819 :
820 : // Add the elements which in this case are 2-node edges that link each pin's nodes vertically.
821 : unsigned int elem_id = elem_sub;
822 8034 : for (unsigned int i = 0; i < _npins; i++)
823 192837 : for (unsigned int iz = 0; iz < _n_cells; iz++)
824 : {
825 185001 : Elem * elem = mesh_base.add_elem(std::make_unique<Edge2>());
826 185001 : elem->subdomain_id() = _pin_block_id;
827 185001 : elem->set_id(elem_id++);
828 185001 : const int indx1 = (_n_cells + 1) * i + iz + node_sub;
829 185001 : const int indx2 = (_n_cells + 1) * i + (iz + 1) + node_sub;
830 185001 : elem->set_node(0, mesh_base.node_ptr(indx1));
831 185001 : elem->set_node(1, mesh_base.node_ptr(indx2));
832 : }
833 :
834 198 : mesh_base.subdomain_name(_pin_block_id) = "fuel_pins";
835 : }
836 :
837 : std::unique_ptr<MeshBase>
838 198 : SCMTriAssemblyMeshGenerator::generate()
839 : {
840 198 : auto mesh_base = buildMeshBaseObject();
841 :
842 : BoundaryInfo & boundary_info = mesh_base->get_boundary_info();
843 198 : mesh_base->set_spatial_dimension(3);
844 198 : mesh_base->reserve_elem(_n_cells * (_n_channels + _npins));
845 198 : mesh_base->reserve_nodes((_n_cells + 1) * (_n_channels + _npins));
846 198 : _nodes.resize(_n_channels);
847 : // Add the points for the give x,y subchannel positions. The grid is hexagonal.
848 : // The grid along
849 : // z is irregular to account for Pin spacers. Store pointers in the _nodes
850 : // array so we can keep track of which points are in which channels.
851 : unsigned int node_id = 0;
852 16662 : for (unsigned int i = 0; i < _n_channels; i++)
853 : {
854 16464 : _nodes[i].reserve(_n_cells + 1);
855 420894 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
856 : {
857 404430 : _nodes[i].push_back(mesh_base->add_point(
858 808860 : Point(_subchannel_position[i][0], _subchannel_position[i][1], _z_grid[iz]), node_id++));
859 : }
860 : }
861 :
862 : // Add the elements which in this case are 2-node edges that link each
863 : // subchannel's nodes vertically.
864 : unsigned int elem_id = 0;
865 16662 : for (unsigned int i = 0; i < _n_channels; i++)
866 : {
867 404430 : for (unsigned int iz = 0; iz < _n_cells; iz++)
868 : {
869 387966 : Elem * elem = mesh_base->add_elem(std::make_unique<Edge2>());
870 387966 : elem->subdomain_id() = _subchannel_block_id;
871 387966 : elem->set_id(elem_id++);
872 387966 : const int indx1 = (_n_cells + 1) * i + iz;
873 387966 : const int indx2 = (_n_cells + 1) * i + (iz + 1);
874 387966 : elem->set_node(0, mesh_base->node_ptr(indx1));
875 387966 : elem->set_node(1, mesh_base->node_ptr(indx2));
876 :
877 387966 : if (iz == 0)
878 16464 : boundary_info.add_side(elem, 0, 0);
879 387966 : if (iz == _n_cells - 1)
880 16464 : boundary_info.add_side(elem, 1, 1);
881 : }
882 : }
883 198 : boundary_info.sideset_name(0) = "inlet";
884 198 : boundary_info.sideset_name(1) = "outlet";
885 198 : boundary_info.nodeset_name(0) = "inlet";
886 198 : boundary_info.nodeset_name(1) = "outlet";
887 :
888 : // Naming the block
889 198 : mesh_base->subdomain_name(_subchannel_block_id) = "subchannel";
890 198 : buildPinMesh(*mesh_base);
891 :
892 198 : mesh_base->prepare_for_use();
893 :
894 : // move the meta data into TriSubChannelMesh
895 198 : auto & sch_mesh = static_cast<TriSubChannelMesh &>(*_mesh);
896 198 : sch_mesh._unheated_length_entry = _unheated_length_entry;
897 198 : sch_mesh._heated_length = _heated_length;
898 198 : sch_mesh._unheated_length_exit = _unheated_length_exit;
899 198 : sch_mesh._z_grid = _z_grid;
900 198 : sch_mesh._k_grid = _k_grid;
901 198 : sch_mesh._spacer_z = _spacer_z;
902 198 : sch_mesh._spacer_k = _spacer_k;
903 198 : sch_mesh._z_blockage = _z_blockage;
904 198 : sch_mesh._index_blockage = _index_blockage;
905 198 : sch_mesh._reduction_blockage = _reduction_blockage;
906 198 : sch_mesh._kij = _kij;
907 198 : sch_mesh._pitch = _pitch;
908 198 : sch_mesh._pin_diameter = _pin_diameter;
909 198 : sch_mesh._n_cells = _n_cells;
910 198 : sch_mesh._n_rings = _n_rings;
911 198 : sch_mesh._n_channels = _n_channels;
912 198 : sch_mesh._flat_to_flat = _flat_to_flat;
913 198 : sch_mesh._dwire = _dwire;
914 198 : sch_mesh._hwire = _hwire;
915 198 : sch_mesh._duct_to_pin_gap = _duct_to_pin_gap;
916 198 : sch_mesh._nodes = _nodes;
917 198 : sch_mesh._gap_to_chan_map = _gap_to_chan_map;
918 198 : sch_mesh._gap_to_pin_map = _gap_to_pin_map;
919 198 : sch_mesh._chan_to_gap_map = _chan_to_gap_map;
920 198 : sch_mesh._sign_id_crossflow_map = _sign_id_crossflow_map;
921 198 : sch_mesh._gij_map = _gij_map;
922 198 : sch_mesh._subchannel_position = _subchannel_position;
923 198 : sch_mesh._pin_position = _pin_position;
924 198 : sch_mesh._pins_in_rings = _pins_in_rings;
925 198 : sch_mesh._chan_to_pin_map = _chan_to_pin_map;
926 198 : sch_mesh._npins = _npins;
927 198 : sch_mesh._n_gaps = _n_gaps;
928 198 : sch_mesh._subch_type = _subch_type;
929 198 : sch_mesh._gap_type = _gap_type;
930 198 : sch_mesh._gap_pairs_sf = _gap_pairs_sf;
931 198 : sch_mesh._chan_pairs_sf = _chan_pairs_sf;
932 198 : sch_mesh._pin_to_chan_map = _pin_to_chan_map;
933 198 : sch_mesh._pin_nodes = _pin_nodes;
934 198 : sch_mesh._pin_mesh_exist = (_npins > 0);
935 198 : sch_mesh.computeAssemblyHydraulicParameters();
936 :
937 198 : return mesh_base;
938 0 : }
|
