Line data Source code
1 : /********************************************************************/
2 : /* SOFTWARE COPYRIGHT NOTIFICATION */
3 : /* Cardinal */
4 : /* */
5 : /* (c) 2021 UChicago Argonne, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by UChicago Argonne, LLC */
9 : /* Under Contract No. DE-AC02-06CH11357 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* Prepared by Battelle Energy Alliance, LLC */
13 : /* Under Contract No. DE-AC07-05ID14517 */
14 : /* With the U. S. Department of Energy */
15 : /* */
16 : /* See LICENSE for full restrictions */
17 : /********************************************************************/
18 :
19 : #include "HexagonalSubchannelMesh.h"
20 : #include "libmesh/cell_prism6.h"
21 : #include "libmesh/face_tri3.h"
22 : #include "libmesh/unstructured_mesh.h"
23 :
24 : registerMooseObject("CardinalApp", HexagonalSubchannelMesh);
25 :
26 : InputParameters
27 208 : HexagonalSubchannelMesh::validParams()
28 : {
29 208 : InputParameters params = HexagonalSubchannelMeshBase::validParams();
30 416 : params.addRequiredRangeCheckedParam<unsigned int>(
31 : "n_axial", "n_axial > 0", "Number of axial cells");
32 416 : params.addRequiredRangeCheckedParam<Real>("height", "height > 0", "Height of assembly");
33 :
34 624 : params.addRangeCheckedParam<unsigned int>(
35 416 : "theta_res", 6, "theta_res >= 2", "Number of nodes on each pin's arc length with a channel");
36 624 : params.addRangeCheckedParam<unsigned int>(
37 416 : "gap_res", 2, "gap_res >= 2", "Number of nodes on each gap");
38 :
39 416 : params.addParam<SubdomainID>("interior_id", 1, "Block ID to set for the interior channels");
40 416 : params.addParam<SubdomainID>("edge_id", 2, "Block ID to set for the edge channels");
41 416 : params.addParam<SubdomainID>("corner_id", 3, "Block ID to set for the corner channels");
42 :
43 416 : params.addParam<bool>("volume_mesh",
44 416 : true,
45 : "Whether to generate a volume mesh (true) "
46 : "or just the surfaces between axial layers in the domain (false)");
47 :
48 208 : params.addClassDescription("Mesh respecting subchannel boundaries for a triangular lattice");
49 208 : return params;
50 0 : }
51 :
52 104 : HexagonalSubchannelMesh::HexagonalSubchannelMesh(const InputParameters & parameters)
53 : : HexagonalSubchannelMeshBase(parameters),
54 104 : _theta_res(getParam<unsigned int>("theta_res")),
55 208 : _gap_res(getParam<unsigned int>("gap_res")),
56 208 : _n_axial(getParam<unsigned int>("n_axial")),
57 208 : _height(getParam<Real>("height")),
58 208 : _interior_id(getParam<SubdomainID>("interior_id")),
59 208 : _edge_id(getParam<SubdomainID>("edge_id")),
60 208 : _corner_id(getParam<SubdomainID>("corner_id")),
61 312 : _volume_mesh(getParam<bool>("volume_mesh"))
62 : {
63 104 : }
64 :
65 : std::unique_ptr<MooseMesh>
66 0 : HexagonalSubchannelMesh::safeClone() const
67 : {
68 0 : return _app.getFactory().copyConstruct(*this);
69 : }
70 :
71 : void
72 104 : HexagonalSubchannelMesh::buildMesh()
73 : {
74 104 : MeshBase & mesh = getMesh();
75 104 : mesh.clear();
76 :
77 104 : _elems_per_interior = 3 * (_theta_res - 1) + 3 * (_gap_res - 1);
78 104 : _elems_per_edge = 2 * (_theta_res - 1) + 4 * (_gap_res - 1);
79 104 : _elems_per_corner = (_theta_res - 1) + 4 * (_gap_res - 1);
80 :
81 104 : _elem_id_counter = 0;
82 104 : _node_id_counter = 0;
83 :
84 104 : Real dz = _height / _n_axial;
85 :
86 : // Build arrays of points corresponding to the three channel types with a centroid
87 : // at (0, 0, 0). These points are then shifted to form the actual channels
88 104 : getInteriorPoints();
89 104 : getEdgePoints();
90 104 : getCornerPoints();
91 :
92 104 : auto nl = _volume_mesh ? _n_axial : _n_axial + 1;
93 :
94 104 : mesh.set_mesh_dimension(3);
95 104 : mesh.set_spatial_dimension(3);
96 :
97 520 : for (unsigned int i = 0; i < nl; ++i)
98 : {
99 416 : Real zmin = i * dz;
100 416 : Real zmax = (i + 1) * dz;
101 :
102 : // add the interior channels
103 416 : if (_hex_lattice.nInteriorChannels())
104 : {
105 : // Then add the elements for the interior channels
106 5516 : for (unsigned int i = 0; i < _hex_lattice.nInteriorChannels(); ++i)
107 : {
108 : Point centroid =
109 5160 : _hex_lattice.channelCentroid(_hex_lattice.interiorChannelCornerCoordinates(i));
110 7740 : Real rotation = i % 2 == 0 ? M_PI : 0.0;
111 :
112 : std::vector<Point> points;
113 114000 : for (const auto & i : _interior_points)
114 108840 : points.push_back(rotatePoint(i, rotation));
115 :
116 108840 : for (int j = 0; j < _elems_per_interior; ++j)
117 : {
118 103680 : bool last_elem = (j == _elems_per_interior - 1);
119 :
120 : Point pt1 = centroid + points[0];
121 103680 : Point pt2 = centroid + points[j + 1];
122 103680 : Point pt3 = last_elem ? centroid + points[1] : centroid + points[j + 2];
123 :
124 103680 : if (_volume_mesh)
125 27216 : addPrismElem(pt1, pt2, pt3, zmin, zmax, _interior_id);
126 : else
127 76464 : addTriElem(pt1, pt2, pt3, zmin, _interior_id);
128 : }
129 5160 : }
130 : }
131 :
132 416 : if (_hex_lattice.nEdgeChannels())
133 : {
134 356 : Real rotation = 0.0;
135 3500 : for (unsigned int i = 0; i < _hex_lattice.nEdgeChannels(); ++i)
136 : {
137 3144 : if (i >= (_n_rings - 1) && i % (_n_rings - 1) == 0)
138 1780 : rotation += 2 * M_PI / 6.0;
139 :
140 3144 : Point centroid = _hex_lattice.channelCentroid(_hex_lattice.edgeChannelCornerCoordinates(i));
141 :
142 : std::vector<Point> points;
143 57792 : for (const auto & i : _edge_points)
144 54648 : points.push_back(rotatePoint(i, rotation));
145 :
146 54648 : for (int j = 0; j < _elems_per_edge; ++j)
147 : {
148 51504 : bool last_elem = (j == _elems_per_edge - 1);
149 :
150 : Point pt1 = centroid + points[0];
151 51504 : Point pt2 = centroid + points[j + 1];
152 51504 : Point pt3 = last_elem ? centroid + points[1] : centroid + points[j + 2];
153 :
154 51504 : if (_volume_mesh)
155 15264 : addPrismElem(pt1, pt2, pt3, zmin, zmax, _edge_id);
156 : else
157 36240 : addTriElem(pt1, pt2, pt3, zmin, _edge_id);
158 : }
159 3144 : }
160 : }
161 :
162 : // there are always corner channels
163 2912 : for (unsigned int i = 0; i < _hex_lattice.nCornerChannels(); ++i)
164 : {
165 2496 : Point centroid = _hex_lattice.channelCentroid(_hex_lattice.cornerChannelCornerCoordinates(i));
166 :
167 : std::vector<Point> points;
168 33768 : for (const auto & pt : _corner_points)
169 31272 : points.push_back(rotatePoint(pt, i * 2 * M_PI / NUM_SIDES));
170 :
171 28776 : for (int j = 0; j < _elems_per_corner; ++j)
172 : {
173 26280 : bool last_elem = (j == _elems_per_corner - 1);
174 :
175 : Point pt1 = centroid + points[0];
176 26280 : Point pt2 = centroid + points[j + 1];
177 26280 : Point pt3 = last_elem ? centroid + points[1] : centroid + points[j + 2];
178 :
179 26280 : if (_volume_mesh)
180 9072 : addPrismElem(pt1, pt2, pt3, zmin, zmax, _corner_id);
181 : else
182 17208 : addTriElem(pt1, pt2, pt3, zmin, _corner_id);
183 : }
184 2496 : }
185 : }
186 :
187 104 : mesh.prepare_for_use();
188 104 : }
189 :
190 : void
191 129912 : HexagonalSubchannelMesh::addTriElem(
192 : const Point & pt1, const Point & pt2, const Point & pt3, const Real & z, const SubdomainID & id)
193 : {
194 129912 : auto elem = new Tri3;
195 129912 : elem->set_id(_elem_id_counter++);
196 129912 : _mesh->add_elem(elem);
197 :
198 129912 : Point z0(0.0, 0.0, z);
199 :
200 129912 : auto node_ptr0 = _mesh->add_point(pt1 + z0, _node_id_counter++);
201 129912 : auto node_ptr1 = _mesh->add_point(pt2 + z0, _node_id_counter++);
202 129912 : auto node_ptr2 = _mesh->add_point(pt3 + z0, _node_id_counter++);
203 :
204 129912 : elem->set_node(0) = node_ptr0;
205 129912 : elem->set_node(1) = node_ptr1;
206 129912 : elem->set_node(2) = node_ptr2;
207 :
208 129912 : elem->subdomain_id() = id;
209 129912 : }
210 :
211 : void
212 51552 : HexagonalSubchannelMesh::addPrismElem(const Point & pt1,
213 : const Point & pt2,
214 : const Point & pt3,
215 : const Real & zmin,
216 : const Real & zmax,
217 : const SubdomainID & id)
218 : {
219 51552 : auto elem = new Prism6;
220 51552 : elem->set_id(_elem_id_counter++);
221 51552 : _mesh->add_elem(elem);
222 :
223 51552 : Point z0(0.0, 0.0, zmin);
224 51552 : Point z1(0.0, 0.0, zmax);
225 :
226 51552 : auto node_ptr0 = _mesh->add_point(pt1 + z0, _node_id_counter++);
227 51552 : auto node_ptr1 = _mesh->add_point(pt2 + z0, _node_id_counter++);
228 51552 : auto node_ptr2 = _mesh->add_point(pt3 + z0, _node_id_counter++);
229 51552 : auto node_ptr3 = _mesh->add_point(pt1 + z1, _node_id_counter++);
230 51552 : auto node_ptr4 = _mesh->add_point(pt2 + z1, _node_id_counter++);
231 51552 : auto node_ptr5 = _mesh->add_point(pt3 + z1, _node_id_counter++);
232 :
233 51552 : elem->set_node(0) = node_ptr0;
234 51552 : elem->set_node(1) = node_ptr1;
235 51552 : elem->set_node(2) = node_ptr2;
236 51552 : elem->set_node(3) = node_ptr3;
237 51552 : elem->set_node(4) = node_ptr4;
238 51552 : elem->set_node(5) = node_ptr5;
239 :
240 51552 : elem->subdomain_id() = id;
241 51552 : }
242 :
243 : void
244 104 : HexagonalSubchannelMesh::getInteriorPoints()
245 : {
246 : int p = 0;
247 104 : _interior_points.resize(3 * _theta_res + 3 * (_gap_res - 2) + 1);
248 104 : _interior_points[p++] = Point(0.0, 0.0, 0.0);
249 :
250 104 : Real distance_from_centroid = _hex_lattice.triangleHeight(_pin_pitch) * 2.0 / 3.0;
251 : Real total_interior_pin_theta = 2.0 * M_PI / 6.0;
252 104 : Real pin_arc_theta = total_interior_pin_theta / (_theta_res - 1.0);
253 104 : Real gap_dx = (_pin_pitch - _pin_diameter) / (_gap_res - 1.0);
254 104 : Real r = _hex_lattice.pinRadius();
255 :
256 : // center coordinates of the three pins forming the channel corners
257 : Point pin1(0.0, distance_from_centroid, 0.0);
258 104 : Point pin2(-distance_from_centroid * COS30, -distance_from_centroid * SIN30, 0.0);
259 : Point pin3(distance_from_centroid * COS30, -distance_from_centroid * SIN30, 0.0);
260 :
261 : // Add the points on the first pin
262 720 : for (unsigned int i = 0; i < _theta_res; ++i)
263 : {
264 616 : Real phi = total_interior_pin_theta + i * pin_arc_theta;
265 616 : _interior_points[p++] = pin1 + Point(r * std::cos(phi), -r * std::sin(phi), 0.0);
266 : }
267 :
268 : // Add the points on the first gap
269 104 : Point start1 = _interior_points[_theta_res];
270 128 : for (unsigned int i = 0; i < _gap_res - 2; ++i)
271 24 : _interior_points[p++] =
272 24 : start1 + Point(-gap_dx * (i + 1) * SIN30, -gap_dx * (i + 1) * COS30, 0.0);
273 :
274 : // Add the points on the second pin
275 720 : for (unsigned int i = 0; i < _theta_res; ++i)
276 : {
277 616 : Real phi = total_interior_pin_theta - i * pin_arc_theta;
278 616 : _interior_points[p++] = pin2 + Point(r * std::cos(phi), r * std::sin(phi), 0.0);
279 : }
280 :
281 : // Add the points on the second gap
282 104 : Point start2 = _interior_points[2 * _theta_res + (_gap_res - 2)];
283 128 : for (unsigned int i = 0; i < _gap_res - 2; ++i)
284 24 : _interior_points[p++] = start2 + Point(gap_dx * (i + 1), 0.0, 0.0);
285 :
286 : // Add points on the third pin
287 720 : for (unsigned int i = 0; i < _theta_res; ++i)
288 : {
289 616 : Real phi = i * pin_arc_theta;
290 616 : _interior_points[p++] = pin3 + Point(-r * std::cos(phi), r * std::sin(phi), 0.0);
291 : }
292 :
293 : // Add points on the third gap
294 104 : Point start3 = _interior_points[3 * _theta_res + 2 * (_gap_res - 2)];
295 128 : for (unsigned int i = 0; i < _gap_res - 2; ++i)
296 24 : _interior_points[p++] =
297 24 : start3 + Point(-gap_dx * (i + 1) * SIN30, gap_dx * (i + 1) * COS30, 0.0);
298 104 : }
299 :
300 : void
301 104 : HexagonalSubchannelMesh::getEdgePoints()
302 : {
303 : int p = 0;
304 104 : _edge_points.resize(2 * _theta_res + 2 * (_gap_res - 2) + 2 * (_gap_res - 1) + 1);
305 :
306 104 : Real width = _pin_pitch;
307 104 : Real r = _hex_lattice.pinRadius();
308 104 : Real height = _hex_lattice.pinBundleSpacing() + r;
309 104 : Real vertical_gap_arc_length = _hex_lattice.pinBundleSpacing() / (_gap_res - 1.0);
310 104 : Real horizontal_gap_arc_length = (_pin_pitch - _pin_diameter) / (_gap_res - 1.0);
311 : Real total_edge_pin_theta = M_PI / 2.0;
312 104 : Real pin_arc_theta = total_edge_pin_theta / (_theta_res - 1.0);
313 :
314 : // center point is at half the distance from the pin outer surface to the duct
315 104 : _edge_points[p++] = Point(0.0, r + _hex_lattice.pinBundleSpacing() / 2.0 - height / 2.0, 0.0);
316 :
317 : // Add points on the first gap
318 104 : Point start1 = Point(-width / 2.0, height / 2.0, 0.0);
319 232 : for (unsigned int i = 0; i < _gap_res - 1; ++i)
320 128 : _edge_points[p++] = start1 + Point(0.0, -vertical_gap_arc_length * i, 0.0);
321 :
322 : // Add points on the first pin
323 : Point pin1 = Point(-width / 2.0, -height / 2.0, 0.0);
324 720 : for (unsigned int i = 0; i < _theta_res; ++i)
325 : {
326 616 : Real phi = pin_arc_theta * i;
327 616 : _edge_points[p++] = pin1 + Point(r * std::sin(phi), r * std::cos(phi), 0.0);
328 : }
329 :
330 : // Add points on the second gap
331 104 : Point start2 = _edge_points[p - 1];
332 128 : for (unsigned int i = 0; i < _gap_res - 2; ++i)
333 24 : _edge_points[p++] = start2 + Point((i + 1) * horizontal_gap_arc_length, 0.0, 0.0);
334 :
335 : // Add points on the second pin
336 : Point pin2 = Point(width / 2.0, -height / 2.0, 0.0);
337 720 : for (unsigned int i = 0; i < _theta_res; ++i)
338 : {
339 616 : Real phi = i * pin_arc_theta;
340 616 : _edge_points[p++] = pin2 + Point(-r * std::cos(phi), r * std::sin(phi), 0.0);
341 : }
342 :
343 : // Add points on the third gap
344 104 : Point start3 = _edge_points[p - 1];
345 232 : for (unsigned int i = 0; i < _gap_res - 1; ++i)
346 128 : _edge_points[p++] = start3 + Point(0.0, vertical_gap_arc_length * (i + 1), 0.0);
347 :
348 : // Add points on the duct
349 104 : Point start4 = _edge_points[p - 1];
350 128 : for (unsigned int i = 0; i < _gap_res - 2; ++i)
351 24 : _edge_points[p++] = start4 + Point(-1.0 * (i + 1) * _pin_pitch / (_gap_res - 1.0), 0.0, 0.0);
352 104 : }
353 :
354 : void
355 104 : HexagonalSubchannelMesh::getCornerPoints()
356 : {
357 : int p = 0;
358 104 : _corner_points.resize(_theta_res + 4 * (_gap_res - 1) + 1);
359 :
360 : // Start with the first corner channel, then we'll translate the centroid (relative
361 : // to the bundle center) to (0, 0, 0)
362 104 : std::vector<Point> pts = _hex_lattice.cornerChannelCornerCoordinates(0);
363 :
364 104 : Point centroid = _hex_lattice.channelCentroid(pts);
365 520 : for (auto & pt : pts)
366 : pt -= centroid;
367 :
368 : // the first point is at the intersection between two lines at the
369 : // halfway point between the pin surface and the duct
370 104 : Real r = _hex_lattice.pinRadius();
371 104 : Real dy = pts[0](1) + r + _hex_lattice.pinBundleSpacing() / 2.0;
372 104 : _corner_points[p++] = Point(dy * SIN30 / COS30, dy, 0.0);
373 :
374 104 : Real gap_arc_length = _hex_lattice.pinBundleSpacing() / (_gap_res - 1.0);
375 : Real total_corner_pin_theta = 2.0 * M_PI / 6.0;
376 104 : Real pin_arc_theta = total_corner_pin_theta / (_theta_res - 1.0);
377 :
378 : // Add the points on the first gap
379 232 : for (unsigned int i = 0; i < _gap_res - 1; ++i)
380 128 : _corner_points[p++] = pts[3] + Point(0.0, -1.0 * i * gap_arc_length, 0.0);
381 :
382 : // Add the points on the pin
383 720 : for (unsigned int i = 0; i < _theta_res; ++i)
384 : {
385 616 : Real phi = i * pin_arc_theta;
386 616 : _corner_points[p++] = pts[0] + Point(r * std::sin(phi), r * std::cos(phi), 0.0);
387 : }
388 :
389 : // Add the points on the second gap
390 232 : for (unsigned int i = 0; i < _gap_res - 1; ++i)
391 : {
392 128 : Real dx = (i + 1) * gap_arc_length;
393 128 : _corner_points[p++] =
394 128 : _corner_points[_theta_res + _gap_res - 1] + Point(dx * COS30, dx * SIN30, 0.0);
395 : }
396 :
397 : // Add the points on the duct
398 : Point wall = pts[2] - pts[1];
399 232 : for (unsigned int i = 0; i < _gap_res - 1; ++i)
400 128 : _corner_points[p++] = pts[1] + wall * (i + 1) / (_gap_res - 1.0);
401 :
402 : wall = pts[3] - pts[2];
403 128 : for (unsigned int i = 0; i < _gap_res - 2; ++i)
404 24 : _corner_points[p++] = pts[2] + wall * (i + 1) / (_gap_res - 1.0);
405 104 : }
|
