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