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 224 : HexagonalSubchannelMesh::validParams()
28 : {
29 224 : InputParameters params = HexagonalSubchannelMeshBase::validParams();
30 448 : params.addRequiredRangeCheckedParam<unsigned int>(
31 : "n_axial", "n_axial > 0", "Number of axial cells");
32 448 : params.addRequiredRangeCheckedParam<Real>("height", "height > 0", "Height of assembly");
33 :
34 672 : params.addRangeCheckedParam<unsigned int>(
35 448 : "theta_res", 6, "theta_res >= 2", "Number of nodes on each pin's arc length with a channel");
36 672 : params.addRangeCheckedParam<unsigned int>(
37 448 : "gap_res", 2, "gap_res >= 2", "Number of nodes on each gap");
38 :
39 448 : params.addParam<SubdomainID>("interior_id", 1, "Block ID to set for the interior channels");
40 448 : params.addParam<SubdomainID>("edge_id", 2, "Block ID to set for the edge channels");
41 448 : params.addParam<SubdomainID>("corner_id", 3, "Block ID to set for the corner channels");
42 :
43 448 : params.addParam<bool>("volume_mesh",
44 448 : true,
45 : "Whether to generate a volume mesh (true) "
46 : "or just the surfaces between axial layers in the domain (false)");
47 :
48 224 : params.addClassDescription("Mesh respecting subchannel boundaries for a triangular lattice");
49 224 : return params;
50 0 : }
51 :
52 112 : HexagonalSubchannelMesh::HexagonalSubchannelMesh(const InputParameters & parameters)
53 : : HexagonalSubchannelMeshBase(parameters),
54 112 : _theta_res(getParam<unsigned int>("theta_res")),
55 224 : _gap_res(getParam<unsigned int>("gap_res")),
56 224 : _n_axial(getParam<unsigned int>("n_axial")),
57 224 : _height(getParam<Real>("height")),
58 224 : _interior_id(getParam<SubdomainID>("interior_id")),
59 224 : _edge_id(getParam<SubdomainID>("edge_id")),
60 224 : _corner_id(getParam<SubdomainID>("corner_id")),
61 336 : _volume_mesh(getParam<bool>("volume_mesh"))
62 : {
63 112 : }
64 :
65 : std::unique_ptr<MooseMesh>
66 0 : HexagonalSubchannelMesh::safeClone() const
67 : {
68 0 : return _app.getFactory().copyConstruct(*this);
69 : }
70 :
71 : void
72 112 : HexagonalSubchannelMesh::buildMesh()
73 : {
74 112 : MeshBase & mesh = getMesh();
75 112 : mesh.clear();
76 :
77 112 : _elems_per_interior = 3 * (_theta_res - 1) + 3 * (_gap_res - 1);
78 112 : _elems_per_edge = 2 * (_theta_res - 1) + 4 * (_gap_res - 1);
79 112 : _elems_per_corner = (_theta_res - 1) + 4 * (_gap_res - 1);
80 :
81 112 : _elem_id_counter = 0;
82 112 : _node_id_counter = 0;
83 :
84 112 : 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 112 : getInteriorPoints();
89 112 : getEdgePoints();
90 112 : getCornerPoints();
91 :
92 112 : auto nl = _volume_mesh ? _n_axial : _n_axial + 1;
93 :
94 112 : mesh.set_mesh_dimension(3);
95 112 : mesh.set_spatial_dimension(3);
96 :
97 584 : for (unsigned int i = 0; i < nl; ++i)
98 : {
99 472 : Real zmin = i * dz;
100 472 : Real zmax = (i + 1) * dz;
101 :
102 : // add the interior channels
103 472 : if (_hex_lattice.nInteriorChannels())
104 : {
105 : // Then add the elements for the interior channels
106 5908 : for (unsigned int i = 0; i < _hex_lattice.nInteriorChannels(); ++i)
107 : {
108 : Point centroid =
109 5496 : _hex_lattice.channelCentroid(_hex_lattice.interiorChannelCornerCoordinates(i));
110 8244 : Real rotation = i % 2 == 0 ? M_PI : 0.0;
111 :
112 : std::vector<Point> points;
113 122736 : for (const auto & i : _interior_points)
114 117240 : points.push_back(rotatePoint(i, rotation));
115 :
116 117240 : for (int j = 0; j < _elems_per_interior; ++j)
117 : {
118 111744 : bool last_elem = (j == _elems_per_interior - 1);
119 :
120 : Point pt1 = centroid + points[0];
121 111744 : Point pt2 = centroid + points[j + 1];
122 111744 : Point pt3 = last_elem ? centroid + points[1] : centroid + points[j + 2];
123 :
124 111744 : if (_volume_mesh)
125 35280 : addPrismElem(pt1, pt2, pt3, zmin, zmax, _interior_id);
126 : else
127 76464 : addTriElem(pt1, pt2, pt3, zmin, _interior_id);
128 : }
129 : }
130 : }
131 :
132 472 : if (_hex_lattice.nEdgeChannels())
133 : {
134 412 : Real rotation = 0.0;
135 3892 : for (unsigned int i = 0; i < _hex_lattice.nEdgeChannels(); ++i)
136 : {
137 3480 : if (i >= (_n_rings - 1) && i % (_n_rings - 1) == 0)
138 2060 : rotation += 2 * M_PI / 6.0;
139 :
140 3480 : Point centroid = _hex_lattice.channelCentroid(_hex_lattice.edgeChannelCornerCoordinates(i));
141 :
142 : std::vector<Point> points;
143 64512 : for (const auto & i : _edge_points)
144 61032 : points.push_back(rotatePoint(i, rotation));
145 :
146 61032 : for (int j = 0; j < _elems_per_edge; ++j)
147 : {
148 57552 : bool last_elem = (j == _elems_per_edge - 1);
149 :
150 : Point pt1 = centroid + points[0];
151 57552 : Point pt2 = centroid + points[j + 1];
152 57552 : Point pt3 = last_elem ? centroid + points[1] : centroid + points[j + 2];
153 :
154 57552 : if (_volume_mesh)
155 21312 : addPrismElem(pt1, pt2, pt3, zmin, zmax, _edge_id);
156 : else
157 36240 : addTriElem(pt1, pt2, pt3, zmin, _edge_id);
158 : }
159 : }
160 : }
161 :
162 : // there are always corner channels
163 3304 : for (unsigned int i = 0; i < _hex_lattice.nCornerChannels(); ++i)
164 : {
165 2832 : Point centroid = _hex_lattice.channelCentroid(_hex_lattice.cornerChannelCornerCoordinates(i));
166 :
167 : std::vector<Point> points;
168 38472 : for (const auto & pt : _corner_points)
169 35640 : points.push_back(rotatePoint(pt, i * 2 * M_PI / NUM_SIDES));
170 :
171 32808 : for (int j = 0; j < _elems_per_corner; ++j)
172 : {
173 29976 : bool last_elem = (j == _elems_per_corner - 1);
174 :
175 : Point pt1 = centroid + points[0];
176 29976 : Point pt2 = centroid + points[j + 1];
177 29976 : Point pt3 = last_elem ? centroid + points[1] : centroid + points[j + 2];
178 :
179 29976 : if (_volume_mesh)
180 12768 : addPrismElem(pt1, pt2, pt3, zmin, zmax, _corner_id);
181 : else
182 17208 : addTriElem(pt1, pt2, pt3, zmin, _corner_id);
183 : }
184 : }
185 : }
186 :
187 112 : mesh.prepare_for_use();
188 112 : }
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 69360 : 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 69360 : auto elem = new Prism6;
220 69360 : elem->set_id(_elem_id_counter++);
221 69360 : _mesh->add_elem(elem);
222 :
223 69360 : Point z0(0.0, 0.0, zmin);
224 69360 : Point z1(0.0, 0.0, zmax);
225 :
226 69360 : auto node_ptr0 = _mesh->add_point(pt1 + z0, _node_id_counter++);
227 69360 : auto node_ptr1 = _mesh->add_point(pt2 + z0, _node_id_counter++);
228 69360 : auto node_ptr2 = _mesh->add_point(pt3 + z0, _node_id_counter++);
229 69360 : auto node_ptr3 = _mesh->add_point(pt1 + z1, _node_id_counter++);
230 69360 : auto node_ptr4 = _mesh->add_point(pt2 + z1, _node_id_counter++);
231 69360 : auto node_ptr5 = _mesh->add_point(pt3 + z1, _node_id_counter++);
232 :
233 69360 : elem->set_node(0) = node_ptr0;
234 69360 : elem->set_node(1) = node_ptr1;
235 69360 : elem->set_node(2) = node_ptr2;
236 69360 : elem->set_node(3) = node_ptr3;
237 69360 : elem->set_node(4) = node_ptr4;
238 69360 : elem->set_node(5) = node_ptr5;
239 :
240 69360 : elem->subdomain_id() = id;
241 69360 : }
242 :
243 : void
244 112 : HexagonalSubchannelMesh::getInteriorPoints()
245 : {
246 : int p = 0;
247 112 : _interior_points.resize(3 * _theta_res + 3 * (_gap_res - 2) + 1);
248 112 : _interior_points[p++] = Point(0.0, 0.0, 0.0);
249 :
250 112 : 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 112 : Real pin_arc_theta = total_interior_pin_theta / (_theta_res - 1.0);
253 112 : Real gap_dx = (_pin_pitch - _pin_diameter) / (_gap_res - 1.0);
254 112 : 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 112 : 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 792 : for (unsigned int i = 0; i < _theta_res; ++i)
263 : {
264 680 : Real phi = total_interior_pin_theta + i * pin_arc_theta;
265 680 : _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 112 : Point start1 = _interior_points[_theta_res];
270 136 : 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 792 : for (unsigned int i = 0; i < _theta_res; ++i)
276 : {
277 680 : Real phi = total_interior_pin_theta - i * pin_arc_theta;
278 680 : _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 112 : Point start2 = _interior_points[2 * _theta_res + (_gap_res - 2)];
283 136 : 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 792 : for (unsigned int i = 0; i < _theta_res; ++i)
288 : {
289 680 : Real phi = i * pin_arc_theta;
290 680 : _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 112 : Point start3 = _interior_points[3 * _theta_res + 2 * (_gap_res - 2)];
295 136 : 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 112 : }
299 :
300 : void
301 112 : HexagonalSubchannelMesh::getEdgePoints()
302 : {
303 : int p = 0;
304 112 : _edge_points.resize(2 * _theta_res + 2 * (_gap_res - 2) + 2 * (_gap_res - 1) + 1);
305 :
306 112 : Real width = _pin_pitch;
307 112 : Real r = _hex_lattice.pinRadius();
308 112 : Real height = _hex_lattice.pinBundleSpacing() + r;
309 112 : Real vertical_gap_arc_length = _hex_lattice.pinBundleSpacing() / (_gap_res - 1.0);
310 112 : Real horizontal_gap_arc_length = (_pin_pitch - _pin_diameter) / (_gap_res - 1.0);
311 : Real total_edge_pin_theta = M_PI / 2.0;
312 112 : 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 112 : _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 112 : Point start1 = Point(-width / 2.0, height / 2.0, 0.0);
319 248 : for (unsigned int i = 0; i < _gap_res - 1; ++i)
320 136 : _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 792 : for (unsigned int i = 0; i < _theta_res; ++i)
325 : {
326 680 : Real phi = pin_arc_theta * i;
327 680 : _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 112 : Point start2 = _edge_points[p - 1];
332 136 : 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 792 : for (unsigned int i = 0; i < _theta_res; ++i)
338 : {
339 680 : Real phi = i * pin_arc_theta;
340 680 : _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 112 : Point start3 = _edge_points[p - 1];
345 248 : for (unsigned int i = 0; i < _gap_res - 1; ++i)
346 136 : _edge_points[p++] = start3 + Point(0.0, vertical_gap_arc_length * (i + 1), 0.0);
347 :
348 : // Add points on the duct
349 112 : Point start4 = _edge_points[p - 1];
350 136 : 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 112 : }
353 :
354 : void
355 112 : HexagonalSubchannelMesh::getCornerPoints()
356 : {
357 : int p = 0;
358 112 : _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 112 : std::vector<Point> pts = _hex_lattice.cornerChannelCornerCoordinates(0);
363 :
364 112 : Point centroid = _hex_lattice.channelCentroid(pts);
365 560 : 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 112 : Real r = _hex_lattice.pinRadius();
371 112 : Real dy = pts[0](1) + r + _hex_lattice.pinBundleSpacing() / 2.0;
372 112 : _corner_points[p++] = Point(dy * SIN30 / COS30, dy, 0.0);
373 :
374 112 : Real gap_arc_length = _hex_lattice.pinBundleSpacing() / (_gap_res - 1.0);
375 : Real total_corner_pin_theta = 2.0 * M_PI / 6.0;
376 112 : Real pin_arc_theta = total_corner_pin_theta / (_theta_res - 1.0);
377 :
378 : // Add the points on the first gap
379 248 : for (unsigned int i = 0; i < _gap_res - 1; ++i)
380 136 : _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 792 : for (unsigned int i = 0; i < _theta_res; ++i)
384 : {
385 680 : Real phi = i * pin_arc_theta;
386 680 : _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 248 : for (unsigned int i = 0; i < _gap_res - 1; ++i)
391 : {
392 136 : Real dx = (i + 1) * gap_arc_length;
393 136 : _corner_points[p++] =
394 136 : _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 248 : for (unsigned int i = 0; i < _gap_res - 1; ++i)
400 136 : _corner_points[p++] = pts[1] + wall * (i + 1) / (_gap_res - 1.0);
401 :
402 : wall = pts[3] - pts[2];
403 136 : 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 112 : }
|
