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 "TriSubChannelMesh.h"
11 : #include <cmath>
12 : #include "libmesh/node.h"
13 :
14 : registerMooseObject("SubChannelApp", TriSubChannelMesh);
15 :
16 : InputParameters
17 442 : TriSubChannelMesh::validParams()
18 : {
19 442 : InputParameters params = SubChannelMesh::validParams();
20 442 : params.addClassDescription("Creates an subchannel mesh container for a triangular "
21 : "lattice arrangement");
22 442 : return params;
23 0 : }
24 :
25 221 : TriSubChannelMesh::TriSubChannelMesh(const InputParameters & params) : SubChannelMesh(params) {}
26 :
27 0 : TriSubChannelMesh::TriSubChannelMesh(const TriSubChannelMesh & other_mesh)
28 : : SubChannelMesh(other_mesh),
29 0 : _n_rings(other_mesh._n_rings),
30 0 : _n_channels(other_mesh._n_channels),
31 0 : _flat_to_flat(other_mesh._flat_to_flat),
32 0 : _dwire(other_mesh._dwire),
33 0 : _hwire(other_mesh._hwire),
34 0 : _duct_to_pin_gap(other_mesh._duct_to_pin_gap),
35 0 : _nodes(other_mesh._nodes),
36 0 : _gap_to_chan_map(other_mesh._gap_to_chan_map),
37 0 : _gap_to_pin_map(other_mesh._gap_to_pin_map),
38 0 : _chan_to_gap_map(other_mesh._chan_to_gap_map),
39 0 : _sign_id_crossflow_map(other_mesh._sign_id_crossflow_map),
40 0 : _gij_map(other_mesh._gij_map),
41 0 : _pin_position(other_mesh._pin_position),
42 0 : _pins_in_rings(other_mesh._pins_in_rings),
43 0 : _chan_to_pin_map(other_mesh._chan_to_pin_map),
44 0 : _npins(other_mesh._npins),
45 0 : _n_gaps(other_mesh._n_gaps),
46 0 : _subch_type(other_mesh._subch_type),
47 0 : _gap_type(other_mesh._gap_type),
48 0 : _gap_pairs_sf(other_mesh._gap_pairs_sf),
49 0 : _chan_pairs_sf(other_mesh._chan_pairs_sf),
50 0 : _pin_to_chan_map(other_mesh._pin_to_chan_map)
51 : {
52 0 : }
53 :
54 : std::unique_ptr<MooseMesh>
55 0 : TriSubChannelMesh::safeClone() const
56 : {
57 0 : return _app.getFactory().copyConstruct(*this);
58 : }
59 :
60 : unsigned int
61 1116432 : TriSubChannelMesh::getSubchannelIndexFromPoint(const Point & p) const
62 : {
63 : /// Function that returns the subchannel index given a point
64 1116432 : return this->channelIndex(p);
65 : }
66 :
67 : unsigned int
68 1557558 : TriSubChannelMesh::channelIndex(const Point & p) const
69 : {
70 : /// Function that returns the subchannel index given a point
71 : /// Determining a channel index given a point
72 : /// Looping over all subchannels to determine the closest one to the point
73 : /// Special treatment for edge and corner subchannels since deformed elements may lead to wrong transfers
74 :
75 : // Distances to determine the closest subchannel
76 : Real distance0 = 1.0e+8; // Dummy distance to keep updating the closes distance found
77 : Real distance1; // Distance to be updated with the current subchannel distance
78 : unsigned int j = 0; // Index to keep track of the closest point to subchannel
79 :
80 : // Projecting point into hexahedral coordinated to determine if the point belongs to a center
81 : // subchannel
82 1557558 : Real distance_outer_ring = _flat_to_flat / 2 - _duct_to_pin_gap - _pin_diameter / 2;
83 1557558 : Real channel_distance = std::sqrt(std::pow(p(0), 2) + std::pow(p(1), 2));
84 1557558 : Real angle = std::abs(std::atan(p(1) / p(0)));
85 1557558 : Real projection_angle =
86 1557558 : angle - libMesh::pi / 6 - std::trunc(angle / (libMesh::pi / 3)) * (libMesh::pi / 3);
87 1557558 : channel_distance = channel_distance * std::cos(projection_angle);
88 :
89 : // Projecting point on top edge to determine if the point is a corner or edge subchannel by x
90 : // coordinate
91 : Real loc_angle = std::atan(p(1) / p(0));
92 1557558 : if (p(0) <= 0)
93 803440 : loc_angle += libMesh::pi;
94 754118 : else if (p(0) >= 0 && p(1) <= 0)
95 391502 : loc_angle += 2 * libMesh::pi;
96 1557558 : Real rem_ang = std::trunc(loc_angle / (libMesh::pi / 3)) * (libMesh::pi / 3) - libMesh::pi / 3;
97 1557558 : Real x_coord_new = (std::cos(-rem_ang) * p(0) - std::sin(-rem_ang) * p(1));
98 1557558 : Real x_lim = (_n_rings - 1) * _pitch / 2.0;
99 :
100 : // looping over all channels
101 137616210 : for (unsigned int i = 0; i < _n_channels; i++)
102 : {
103 :
104 136058652 : if (_n_rings == 1)
105 : {
106 : Real angle = std::atan(p(1) / p(0));
107 0 : if ((i * libMesh::pi / 6.0 < angle) && (angle <= (i + 1) * libMesh::pi / 6.0))
108 0 : return i;
109 : }
110 : else
111 : {
112 : // Distance from the point to subchannel
113 136058652 : distance1 = std::sqrt(std::pow((p(0) - _subchannel_position[i][0]), 2.0) +
114 136058652 : std::pow((p(1) - _subchannel_position[i][1]), 2.0));
115 :
116 : // If subchannel belongs to center ring
117 136058652 : if (channel_distance < distance_outer_ring)
118 : {
119 96831324 : if ((distance1 < distance0) && (_subch_type[i] == EChannelType::CENTER))
120 : {
121 : j = i;
122 : distance0 = distance1;
123 : } // if
124 : } // if
125 : // If subchannel belongs to outer ring
126 : else
127 : {
128 39227328 : if ((distance1 < distance0) &&
129 30627787 : (_subch_type[i] == EChannelType::EDGE || _subch_type[i] == EChannelType::CORNER))
130 : {
131 3213499 : if (((x_coord_new > x_lim) || (x_coord_new < -x_lim)) &&
132 : _subch_type[i] == EChannelType::CORNER)
133 : {
134 : j = i;
135 : distance0 = distance1;
136 : } // if
137 3013137 : else if (((x_coord_new > x_lim) || (x_coord_new > -x_lim)) &&
138 : _subch_type[i] == EChannelType::EDGE)
139 : {
140 : j = i;
141 : distance0 = distance1;
142 : }
143 : }
144 : }
145 : }
146 : }
147 : return j;
148 : }
149 :
150 : void
151 221 : TriSubChannelMesh::buildMesh()
152 : {
153 221 : }
154 :
155 : void
156 207 : TriSubChannelMesh::computeAssemblyHydraulicParameters()
157 : {
158 207 : _assembly_flow_area = 0.0;
159 207 : _assembly_wetted_perimeter = 0.0;
160 207 : _assembly_hydraulic_diameter = 0.0;
161 :
162 207 : const Real z = _z_grid.empty() ? 0.0 : _z_grid.front();
163 :
164 17505 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
165 : {
166 17298 : _assembly_flow_area += getSubchannelFlowArea(i_ch, z);
167 17298 : _assembly_wetted_perimeter += getSubchannelWettedPerimeter(i_ch);
168 : }
169 :
170 207 : if (_assembly_wetted_perimeter == 0.0)
171 0 : mooseError(name(), ": Assembly wetted perimeter is zero; cannot compute hydraulic diameter.");
172 :
173 207 : _assembly_hydraulic_diameter = 4.0 * _assembly_flow_area / _assembly_wetted_perimeter;
174 207 : }
175 :
176 : Real
177 484818 : TriSubChannelMesh::getSubchannelFlowArea(unsigned int i_chan, Real z) const
178 : {
179 : Real standard_area = 0.0;
180 : Real rod_area = 0.0;
181 : Real wire_area = 0.0;
182 :
183 : const Real theta =
184 484818 : std::acos(_hwire / std::sqrt(std::pow(_hwire, 2) +
185 484818 : std::pow(libMesh::pi * (_pin_diameter + _dwire), 2)));
186 484818 : const auto subch_type = getSubchannelType(i_chan);
187 484818 : if (subch_type == EChannelType::CENTER)
188 : {
189 336624 : standard_area = std::pow(_pitch, 2) * std::sqrt(3.0) / 4.0;
190 336624 : rod_area = libMesh::pi * std::pow(_pin_diameter, 2) / 8.0;
191 336624 : wire_area = libMesh::pi * std::pow(_dwire, 2) / 8.0 / std::cos(theta);
192 : }
193 148194 : else if (subch_type == EChannelType::EDGE)
194 : {
195 108072 : standard_area = _pitch * (_pin_diameter / 2.0 + _duct_to_pin_gap);
196 108072 : rod_area = libMesh::pi * std::pow(_pin_diameter, 2) / 8.0;
197 108072 : wire_area = libMesh::pi * std::pow(_dwire, 2) / 8.0 / std::cos(theta);
198 : }
199 : else
200 : {
201 40122 : standard_area = 1.0 / std::sqrt(3.0) * std::pow(_pin_diameter / 2.0 + _duct_to_pin_gap, 2);
202 40122 : rod_area = libMesh::pi * std::pow(_pin_diameter, 2) / 24.0;
203 40122 : wire_area = libMesh::pi * std::pow(_dwire, 2) / 24.0 / std::cos(theta);
204 : }
205 :
206 484818 : Real flow_area = standard_area - rod_area - wire_area;
207 :
208 : unsigned int blockage_index = 0;
209 1600350 : for (const auto & i_blockage : _index_blockage)
210 : {
211 1115532 : if (i_chan == i_blockage && z >= _z_blockage.front() && z <= _z_blockage.back())
212 638 : flow_area *= _reduction_blockage[blockage_index];
213 1115532 : blockage_index++;
214 : }
215 :
216 484818 : return flow_area;
217 : }
218 :
219 : Real
220 484818 : TriSubChannelMesh::getSubchannelWettedPerimeter(unsigned int i_chan) const
221 : {
222 : const Real theta =
223 484818 : std::acos(_hwire / std::sqrt(std::pow(_hwire, 2) +
224 484818 : std::pow(libMesh::pi * (_pin_diameter + _dwire), 2)));
225 484818 : const Real rod_circumference = libMesh::pi * _pin_diameter;
226 484818 : const Real wire_circumference = libMesh::pi * _dwire;
227 484818 : const auto subch_type = getSubchannelType(i_chan);
228 :
229 484818 : if (subch_type == EChannelType::CENTER)
230 336624 : return 0.5 * rod_circumference + 0.5 * wire_circumference / std::cos(theta);
231 148194 : else if (subch_type == EChannelType::EDGE)
232 108072 : return 0.5 * rod_circumference + 0.5 * wire_circumference / std::cos(theta) + _pitch;
233 : else
234 40122 : return (rod_circumference + wire_circumference / std::cos(theta)) / 6.0 +
235 40122 : 2.0 / std::sqrt(3.0) * (_pin_diameter / 2.0 + _duct_to_pin_gap);
236 : }
237 :
238 : unsigned int
239 92613 : TriSubChannelMesh::getPinIndexFromPoint(const Point & p) const
240 : {
241 : /// Function that returns the pin number given a point
242 :
243 92613 : return this->pinIndex(p);
244 : }
245 :
246 : unsigned int
247 92613 : TriSubChannelMesh::pinIndex(const Point & p) const
248 : {
249 : /// Function that returns the pin number given a point
250 :
251 : // Define the current ring
252 : Real distance_rod;
253 : Real d0 = 1e5;
254 : unsigned int current_rod = 0;
255 :
256 : std::vector<Point> positions;
257 : Point center(0, 0);
258 92613 : this->pinPositions(positions, _n_rings, _pitch, center);
259 7136286 : for (unsigned int i = 0; i < _npins; i++)
260 : {
261 7043673 : Real x_dist = positions[i](0) - p(0);
262 7043673 : Real y_dist = positions[i](1) - p(1);
263 7043673 : distance_rod = std::sqrt(std::pow(x_dist, 2) + std::pow(y_dist, 2));
264 7043673 : if (distance_rod < d0)
265 : {
266 : d0 = distance_rod;
267 : current_rod = i;
268 : }
269 : }
270 :
271 92613 : return current_rod;
272 92613 : }
273 :
274 : void
275 92953 : TriSubChannelMesh::pinPositions(std::vector<Point> & positions,
276 : unsigned int nrings,
277 : Real pitch,
278 : Point center)
279 : {
280 : /// Defining parameters
281 : // distance: it is the distance to the next Pin
282 : //
283 : Real theta = 0.0;
284 : Real dtheta = 0.0;
285 : Real distance = 0.0;
286 : Real theta1 = 0.0;
287 : Real theta_corrected = 0.0;
288 : Real pi = libMesh::pi;
289 : unsigned int k = 0;
290 92953 : positions.emplace_back(0.0, 0.0);
291 494605 : for (unsigned int i = 1; i < nrings; i++)
292 : {
293 401652 : dtheta = 2.0 * pi / (i * 6);
294 : theta = 0.0;
295 7367118 : for (unsigned int j = 0; j < i * 6; j++)
296 : {
297 : k = k + 1;
298 6965466 : theta1 = fmod(theta + 1.0e-10, pi / 3.0);
299 6965466 : distance = std::sqrt((pow(i * pitch, 2) + pow(theta1 / dtheta * pitch, 2) -
300 6965466 : 2.0 * i * pitch * (theta1 / dtheta * pitch) * std::cos(pi / 3.0)));
301 6965466 : double argument = 1.0 / (i * pitch) / distance / 2.0 *
302 6965466 : (std::pow(i * pitch, 2.0) + std::pow(distance, 2.0) -
303 6965466 : std::pow(theta1 / dtheta * pitch, 2.0));
304 : // Check if the argument to std::acos() is within the valid range [-1, 1]
305 6965466 : if (argument >= -1.0 && argument <= 1.0)
306 : {
307 6754137 : theta_corrected = std::acos(argument);
308 : }
309 : else if (argument > 1.0)
310 : {
311 : theta_corrected = 0.0;
312 : }
313 : else
314 : {
315 : theta_corrected = pi;
316 : }
317 6965466 : if (theta1 < 1.0e-6)
318 : {
319 : theta_corrected = theta;
320 : }
321 : else
322 : {
323 4555554 : if (theta > pi / 3.0 && theta <= 2.0 / 3.0 * pi)
324 759259 : theta_corrected = theta_corrected + pi / 3.0;
325 3037036 : else if (theta > 2.0 / 3.0 * pi && theta <= pi)
326 759259 : theta_corrected = theta_corrected + 2.0 / 3.0 * pi;
327 2277777 : else if (theta > pi && theta <= 4.0 / 3.0 * pi)
328 759259 : theta_corrected = theta_corrected + pi;
329 1518518 : else if (theta > 4.0 / 3.0 * pi && theta <= 5.0 / 3.0 * pi)
330 759259 : theta_corrected = theta_corrected + 4.0 / 3.0 * pi;
331 759259 : else if (theta > 5.0 / 3.0 * pi && theta <= 2.0 * pi)
332 759259 : theta_corrected = theta_corrected + 5.0 / 3.0 * pi;
333 : }
334 13930932 : positions.emplace_back(center(0) + distance * std::cos(theta_corrected),
335 6965466 : center(1) + distance * std::sin(theta_corrected));
336 6965466 : theta = theta + dtheta;
337 : } // j
338 : } // i
339 92953 : }
|
