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 "TriSubChannel1PhaseProblem.h"
11 : #include "AuxiliarySystem.h"
12 : #include "TriSubChannelMesh.h"
13 : #include "SCM.h"
14 : #include <limits> // for std::numeric_limits
15 : #include <cmath> // for std::isnan
16 :
17 : registerMooseObject("SubChannelApp", TriSubChannel1PhaseProblem);
18 :
19 : InputParameters
20 316 : TriSubChannel1PhaseProblem::validParams()
21 : {
22 316 : InputParameters params = SubChannel1PhaseProblem::validParams();
23 316 : params.addClassDescription("Solver class for subchannels in a triangular lattice assembly and "
24 : "bare/wire-wrapped fuel pins");
25 316 : return params;
26 0 : }
27 :
28 158 : TriSubChannel1PhaseProblem::TriSubChannel1PhaseProblem(const InputParameters & params)
29 : : SubChannel1PhaseProblem(params),
30 158 : _tri_sch_mesh(SCM::getMesh<TriSubChannelMesh>(_subchannel_mesh))
31 : {
32 : // Initializing heat conduction system
33 158 : LibmeshPetscCall(createPetscMatrix(
34 : _hc_axial_heat_conduction_mat, _block_size * _n_channels, _block_size * _n_channels));
35 158 : LibmeshPetscCall(createPetscVector(_hc_axial_heat_conduction_rhs, _block_size * _n_channels));
36 158 : LibmeshPetscCall(createPetscMatrix(
37 : _hc_radial_heat_conduction_mat, _block_size * _n_channels, _block_size * _n_channels));
38 158 : LibmeshPetscCall(createPetscVector(_hc_radial_heat_conduction_rhs, _block_size * _n_channels));
39 158 : LibmeshPetscCall(createPetscMatrix(
40 : _hc_sweep_enthalpy_mat, _block_size * _n_channels, _block_size * _n_channels));
41 158 : LibmeshPetscCall(createPetscVector(_hc_sweep_enthalpy_rhs, _block_size * _n_channels));
42 158 : }
43 :
44 158 : TriSubChannel1PhaseProblem::~TriSubChannel1PhaseProblem()
45 : {
46 158 : PetscErrorCode ierr = cleanUp();
47 158 : if (ierr)
48 0 : mooseError(name(), ": Error in memory cleanup");
49 158 : }
50 :
51 : PetscErrorCode
52 158 : TriSubChannel1PhaseProblem::cleanUp()
53 : {
54 : PetscFunctionBegin;
55 : // Clean up heat conduction system
56 158 : LibmeshPetscCall(MatDestroy(&_hc_axial_heat_conduction_mat));
57 158 : LibmeshPetscCall(VecDestroy(&_hc_axial_heat_conduction_rhs));
58 158 : LibmeshPetscCall(MatDestroy(&_hc_radial_heat_conduction_mat));
59 158 : LibmeshPetscCall(VecDestroy(&_hc_radial_heat_conduction_rhs));
60 158 : LibmeshPetscCall(MatDestroy(&_hc_sweep_enthalpy_mat));
61 158 : LibmeshPetscCall(VecDestroy(&_hc_sweep_enthalpy_rhs));
62 158 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
63 : }
64 :
65 : void
66 153 : TriSubChannel1PhaseProblem::initializeSolution()
67 : {
68 153 : if (_deformation)
69 : {
70 : // update surface area, wetted perimeter based on: Dpin, displacement
71 : Real standard_area, wire_area, additional_area, wetted_perimeter, displaced_area;
72 9 : auto flat_to_flat = _tri_sch_mesh.getFlatToFlat();
73 9 : auto n_rings = _tri_sch_mesh.getNumOfRings();
74 9 : auto pitch = _subchannel_mesh.getPitch();
75 9 : auto pin_diameter = _subchannel_mesh.getPinDiameter();
76 9 : auto wire_diameter = _tri_sch_mesh.getWireDiameter();
77 9 : auto wire_lead_length = _tri_sch_mesh.getWireLeadLength();
78 9 : auto gap = _tri_sch_mesh.getDuctToPinGap();
79 9 : auto z_blockage = _subchannel_mesh.getZBlockage();
80 9 : auto index_blockage = _subchannel_mesh.getIndexBlockage();
81 9 : auto reduction_blockage = _subchannel_mesh.getReductionBlockage();
82 9 : auto theta = std::acos(wire_lead_length /
83 9 : std::sqrt(std::pow(wire_lead_length, 2) +
84 9 : std::pow(libMesh::pi * (pin_diameter + wire_diameter), 2)));
85 318 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
86 : {
87 38865 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
88 : {
89 38556 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
90 38556 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
91 38556 : auto Z = _z_grid[iz];
92 : Real rod_area = 0.0;
93 : Real rod_perimeter = 0.0;
94 143208 : for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
95 : {
96 104652 : auto * pin_node = _subchannel_mesh.getPinNode(i_pin, iz);
97 104652 : if (subch_type == EChannelType::CENTER || subch_type == EChannelType::CORNER)
98 : {
99 89964 : rod_area +=
100 89964 : (1.0 / 6.0) * 0.25 * M_PI * (*_Dpin_soln)(pin_node) * (*_Dpin_soln)(pin_node);
101 89964 : rod_perimeter += (1.0 / 6.0) * M_PI * (*_Dpin_soln)(pin_node);
102 : }
103 : else
104 : {
105 14688 : rod_area +=
106 14688 : (1.0 / 4.0) * 0.25 * M_PI * (*_Dpin_soln)(pin_node) * (*_Dpin_soln)(pin_node);
107 14688 : rod_perimeter += (1.0 / 4.0) * M_PI * (*_Dpin_soln)(pin_node);
108 : }
109 : }
110 :
111 38556 : if (subch_type == EChannelType::CENTER)
112 : {
113 29376 : standard_area = std::pow(pitch, 2.0) * std::sqrt(3.0) / 4.0;
114 : additional_area = 0.0;
115 : displaced_area = 0.0;
116 29376 : wire_area = libMesh::pi * std::pow(wire_diameter, 2.0) / 8.0 / std::cos(theta);
117 29376 : wetted_perimeter = rod_perimeter + 0.5 * libMesh::pi * wire_diameter / std::cos(theta);
118 : }
119 9180 : else if (subch_type == EChannelType::EDGE)
120 : {
121 7344 : standard_area = pitch * (pin_diameter / 2.0 + gap);
122 : additional_area = 0.0;
123 7344 : displaced_area = (*_displacement_soln)(node)*pitch;
124 7344 : wire_area = libMesh::pi * std::pow(wire_diameter, 2.0) / 8.0 / std::cos(theta);
125 7344 : wetted_perimeter =
126 7344 : rod_perimeter + 0.5 * libMesh::pi * wire_diameter / std::cos(theta) + pitch;
127 : }
128 : else
129 : {
130 1836 : standard_area = 1.0 / std::sqrt(3.0) * std::pow((pin_diameter / 2.0 + gap), 2.0);
131 : additional_area = 0.0;
132 3672 : displaced_area = 1.0 / std::sqrt(3.0) *
133 1836 : (pin_diameter + 2.0 * gap + (*_displacement_soln)(node)) *
134 1836 : (*_displacement_soln)(node);
135 1836 : wire_area = libMesh::pi / 24.0 * std::pow(wire_diameter, 2.0) / std::cos(theta);
136 1836 : wetted_perimeter =
137 1836 : rod_perimeter + libMesh::pi * wire_diameter / std::cos(theta) / 6.0 +
138 1836 : 2.0 / std::sqrt(3.0) * (pin_diameter / 2.0 + gap + (*_displacement_soln)(node));
139 : }
140 :
141 : // Calculate subchannel area
142 38556 : auto subchannel_area =
143 38556 : standard_area + additional_area + displaced_area - rod_area - wire_area;
144 :
145 : // Correct subchannel area and wetted perimeter in case of overlapping pins
146 : auto overlapping_pin_area = 0.0;
147 : auto overlapping_wetted_perimeter = 0.0;
148 152388 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
149 : {
150 113832 : auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
151 : auto pin_1 = gap_pins.first;
152 : auto pin_2 = gap_pins.second;
153 113832 : auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
154 113832 : auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
155 113832 : auto Diameter1 = (*_Dpin_soln)(pin_node_1);
156 113832 : auto Radius1 = Diameter1 / 2.0;
157 113832 : auto Diameter2 = (*_Dpin_soln)(pin_node_2);
158 113832 : auto Radius2 = Diameter2 / 2.0;
159 113832 : auto pitch = _subchannel_mesh.getPitch();
160 :
161 113832 : if (pitch < (Radius1 + Radius2)) // overlapping pins
162 : {
163 0 : mooseWarning(" The gap of index : '", i_gap, " at axial cell ", iz, " ' is blocked.");
164 0 : auto cos1 =
165 0 : (pitch * pitch + Radius1 * Radius1 - Radius2 * Radius2) / (2 * pitch * Radius1);
166 0 : auto cos2 =
167 0 : (pitch * pitch + Radius2 * Radius2 - Radius1 * Radius1) / (2 * pitch * Radius2);
168 0 : auto angle1 = 2.0 * acos(cos1);
169 0 : auto angle2 = 2.0 * acos(cos2);
170 : // half of the intersecting arc-length
171 0 : overlapping_wetted_perimeter += 0.5 * angle1 * Radius1 + 0.5 * angle2 * Radius2;
172 : // Half of the overlapping area
173 0 : overlapping_pin_area +=
174 0 : 0.5 * Radius1 * Radius1 * acos(cos1) + 0.5 * Radius2 * Radius2 * acos(cos2) -
175 0 : 0.25 * sqrt((-pitch + Radius1 + Radius2) * (pitch + Radius1 - Radius2) *
176 0 : (pitch - Radius1 + Radius2) * (pitch + Radius1 + Radius2));
177 : }
178 : }
179 38556 : subchannel_area += overlapping_pin_area; // correct surface area
180 38556 : wetted_perimeter += -overlapping_wetted_perimeter; // correct wetted perimeter
181 :
182 : // Apply area reduction on subchannels affected by blockage
183 : auto index = 0;
184 77112 : for (const auto & i_blockage : index_blockage)
185 : {
186 38556 : if (i_ch == i_blockage && (Z >= z_blockage.front() && Z <= z_blockage.back()))
187 : {
188 6 : subchannel_area *= reduction_blockage[index];
189 : }
190 38556 : index++;
191 : }
192 38556 : _S_flow_soln->set(node, subchannel_area);
193 38556 : _w_perim_soln->set(node, wetted_perimeter);
194 : }
195 : }
196 : // update map of gap between pins (gij) based on: Dpin, displacement
197 318 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
198 : {
199 57225 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
200 : {
201 56916 : auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
202 : auto pin_1 = gap_pins.first;
203 : auto pin_2 = gap_pins.second;
204 56916 : auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
205 56916 : auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
206 :
207 56916 : if (pin_1 == pin_2) // Corner or edge gap
208 : {
209 : auto displacement = 0.0;
210 : auto counter = 0.0;
211 55080 : for (auto i_ch : _subchannel_mesh.getPinChannels(pin_1))
212 : {
213 45900 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
214 45900 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
215 45900 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
216 : {
217 22032 : displacement += (*_displacement_soln)(node);
218 22032 : counter += 1.0;
219 : }
220 : }
221 9180 : displacement = displacement / counter;
222 9180 : _tri_sch_mesh._gij_map[iz][i_gap] =
223 9180 : 0.5 * (flat_to_flat - (n_rings - 1) * pitch * std::sqrt(3.0) -
224 9180 : (*_Dpin_soln)(pin_node_1)) +
225 : displacement;
226 : }
227 : else // center gap
228 : {
229 47736 : _tri_sch_mesh._gij_map[iz][i_gap] =
230 47736 : pitch - (*_Dpin_soln)(pin_node_1) / 2.0 - (*_Dpin_soln)(pin_node_2) / 2.0;
231 : }
232 : // if pins come in contact, the gap is zero
233 56916 : if (_tri_sch_mesh._gij_map[iz][i_gap] <= 0.0)
234 0 : _tri_sch_mesh._gij_map[iz][i_gap] = 0.0;
235 : }
236 : }
237 9 : }
238 :
239 3261 : for (unsigned int iz = 1; iz < _n_cells + 1; iz++)
240 : {
241 181164 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
242 : {
243 178056 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
244 178056 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
245 178056 : _mdot_soln->set(node_out, (*_mdot_soln)(node_in));
246 : }
247 : }
248 :
249 : // We must do a global assembly to make sure data is parallel consistent before we do things
250 : // like compute L2 norms
251 153 : _aux->solution().close();
252 153 : }
253 :
254 : Real
255 6409224 : TriSubChannel1PhaseProblem::computeFrictionFactor(FrictionStruct friction_args)
256 : {
257 : // The upgraded Cheng and Todreas correlation for pressure drop in hexagonal wire-wrapped rod
258 : // bundles
259 6409224 : auto Re = friction_args.Re;
260 6409224 : auto i_ch = friction_args.i_ch;
261 6409224 : auto S = friction_args.S;
262 6409224 : auto w_perim = friction_args.w_perim;
263 6409224 : auto Dh_i = 4.0 * S / w_perim;
264 : Real aL, b1L, b2L, cL;
265 : Real aT, b1T, b2T, cT;
266 6409224 : const Real & pitch = _subchannel_mesh.getPitch();
267 6409224 : const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
268 6409224 : const Real & wire_lead_length = _tri_sch_mesh.getWireLeadLength();
269 6409224 : const Real & wire_diameter = _tri_sch_mesh.getWireDiameter();
270 6409224 : auto p_over_d = pitch / pin_diameter;
271 6409224 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
272 : // This gap is a constant value for the whole assembly. Might want to make it
273 : // subchannel specific in the future if we have duct deformation.
274 6409224 : auto gap = _tri_sch_mesh.getDuctToPinGap();
275 6409224 : auto w_over_d = (pin_diameter + gap) / pin_diameter;
276 6409224 : auto ReL = std::pow(10, (p_over_d - 1)) * 320.0;
277 6409224 : auto ReT = std::pow(10, 0.7 * (p_over_d - 1)) * 1.0E+4;
278 6409224 : auto psi = std::log(Re / ReL) / std::log(ReT / ReL);
279 : const Real lambda = 7.0;
280 6409224 : auto theta = std::acos(wire_lead_length /
281 6409224 : std::sqrt(std::pow(wire_lead_length, 2) +
282 6409224 : std::pow(libMesh::pi * (pin_diameter + wire_diameter), 2)));
283 6409224 : auto wd_t = (19.56 - 98.71 * (wire_diameter / pin_diameter) +
284 6409224 : 303.47 * std::pow((wire_diameter / pin_diameter), 2.0)) *
285 6409224 : std::pow((wire_lead_length / pin_diameter), -0.541);
286 6409224 : auto wd_l = 1.4 * wd_t;
287 6409224 : auto ws_t = -11.0 * std::log(wire_lead_length / pin_diameter) + 19.0;
288 : auto ws_l = ws_t;
289 : Real pw_p = 0.0;
290 : Real ar = 0.0;
291 : Real a_p = 0.0;
292 :
293 : // Find the coefficients of bare Pin bundle friction factor
294 : // correlations for turbulent and laminar flow regimes. Todreas & Kazimi, Nuclear Systems
295 : // second edition, Volume 1, Chapter 9.6
296 6409224 : if (subch_type == EChannelType::CENTER)
297 : {
298 3825648 : if (p_over_d < 1.1)
299 : {
300 : aL = 26.0;
301 : b1L = 888.2;
302 : b2L = -3334.0;
303 : aT = 0.09378;
304 : b1T = 1.398;
305 : b2T = -8.664;
306 : }
307 : else
308 : {
309 : aL = 62.97;
310 : b1L = 216.9;
311 : b2L = -190.2;
312 : aT = 0.1458;
313 : b1T = 0.03632;
314 : b2T = -0.03333;
315 : }
316 : // laminar flow friction factor for bare Pin bundle - Center subchannel
317 3825648 : cL = aL + b1L * (p_over_d - 1) + b2L * std::pow((p_over_d - 1), 2.0);
318 : // turbulent flow friction factor for bare Pin bundle - Center subchannel
319 3825648 : cT = aT + b1T * (p_over_d - 1) + b2T * std::pow((p_over_d - 1), 2.0);
320 : }
321 2583576 : else if (subch_type == EChannelType::EDGE)
322 : {
323 1750464 : if (w_over_d < 1.1)
324 : {
325 : aL = 26.18;
326 : b1L = 554.5;
327 : b2L = -1480.0;
328 : aT = 0.09377;
329 : b1T = 0.8732;
330 : b2T = -3.341;
331 : }
332 : else
333 : {
334 : aL = 44.4;
335 : b1L = 256.7;
336 : b2L = -267.6;
337 : aT = 0.1430;
338 : b1T = 0.04199;
339 : b2T = -0.04428;
340 : }
341 : // laminar flow friction factor for bare Pin bundle - Edge subchannel
342 1750464 : cL = aL + b1L * (w_over_d - 1) + b2L * std::pow((w_over_d - 1), 2.0);
343 : // turbulent flow friction factor for bare Pin bundle - Edge subchannel
344 1750464 : cT = aT + b1T * (w_over_d - 1) + b2T * std::pow((w_over_d - 1), 2.0);
345 : }
346 : else
347 : {
348 833112 : if (w_over_d < 1.1)
349 : {
350 : aL = 26.98;
351 : b1L = 1636.0;
352 : b2L = -10050.0;
353 : aT = 0.1004;
354 : b1T = 1.625;
355 : b2T = -11.85;
356 : }
357 : else
358 : {
359 : aL = 87.26;
360 : b1L = 38.59;
361 : b2L = -55.12;
362 : aT = 0.1499;
363 : b1T = 0.006706;
364 : b2T = -0.009567;
365 : }
366 : // laminar flow friction factor for bare Pin bundle - Corner subchannel
367 833112 : cL = aL + b1L * (w_over_d - 1) + b2L * std::pow((w_over_d - 1), 2.0);
368 : // turbulent flow friction factor for bare Pin bundle - Corner subchannel
369 833112 : cT = aT + b1T * (w_over_d - 1) + b2T * std::pow((w_over_d - 1), 2.0);
370 : }
371 :
372 : // Find the coefficients of wire-wrapped Pin bundle friction factor
373 : // correlations for turbulent and laminar flow regimes. Todreas & Kazimi, Nuclear Systems
374 : // Volume 1 Chapter 9-6 also Chen and Todreas (2018).
375 6409224 : if ((wire_diameter != 0.0) && (wire_lead_length != 0.0))
376 : {
377 6128784 : if (subch_type == EChannelType::CENTER)
378 : {
379 : // wetted perimeter for center subchannel and bare Pin bundle
380 3653568 : pw_p = libMesh::pi * pin_diameter / 2.0;
381 : // wire projected area - center subchannel wire-wrapped bundle
382 3653568 : ar = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
383 : // bare Pin bundle center subchannel flow area (normal area + wire area)
384 3653568 : a_p = S + libMesh::pi * std::pow(wire_diameter, 2.0) / 8.0 / std::cos(theta);
385 : // turbulent friction factor equation constant - Center subchannel
386 3653568 : cT *= (pw_p / w_perim);
387 3653568 : cT += wd_t * (3.0 * ar / a_p) * (Dh_i / wire_lead_length) *
388 3653568 : std::pow((Dh_i / wire_diameter), 0.18);
389 : // laminar friction factor equation constant - Center subchannel
390 3653568 : cL *= (pw_p / w_perim);
391 3653568 : cL += wd_l * (3.0 * ar / a_p) * (Dh_i / wire_lead_length) * (Dh_i / wire_diameter);
392 : }
393 2475216 : else if (subch_type == EChannelType::EDGE)
394 : {
395 : // wire projected area - edge subchannel wire-wrapped bundle
396 1677024 : ar = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 4.0;
397 : // bare Pin bundle edge subchannel flow area (normal area + wire area)
398 1677024 : a_p = S + libMesh::pi * std::pow(wire_diameter, 2.0) / 8.0 / std::cos(theta);
399 : // turbulent friction factor equation constant - Edge subchannel
400 1677024 : cT *= std::pow((1 + ws_t * (ar / a_p) * std::pow(std::tan(theta), 2.0)), 1.41);
401 : // laminar friction factor equation constant - Edge subchannel
402 1677024 : cL *= (1 + ws_l * (ar / a_p) * std::pow(std::tan(theta), 2.0));
403 : }
404 : else
405 : {
406 : // wire projected area - corner subchannel wire-wrapped bundle
407 798192 : ar = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
408 : // bare Pin bundle corner subchannel flow area (normal area + wire area)
409 798192 : a_p = S + libMesh::pi * std::pow(wire_diameter, 2.0) / 24.0 / std::cos(theta);
410 : // turbulent friction factor equation constant - Corner subchannel
411 798192 : cT *= std::pow((1 + ws_t * (ar / a_p) * std::pow(std::tan(theta), 2.0)), 1.41);
412 : // laminar friction factor equation constant - Corner subchannel
413 798192 : cL *= (1 + ws_l * (ar / a_p) * std::pow(std::tan(theta), 2.0));
414 : }
415 : }
416 :
417 : // laminar friction factor
418 6409224 : auto fL = cL / Re;
419 : // turbulent friction factor
420 6409224 : auto fT = cT / std::pow(Re, 0.18);
421 :
422 6409224 : if (Re < ReL)
423 : {
424 : // laminar flow
425 : return fL;
426 : }
427 6409224 : else if (Re > ReT)
428 : {
429 : // turbulent flow
430 : return fT;
431 : }
432 : else
433 : {
434 : // transient flow: psi definition uses a Bulk ReT/ReL number, same for all channels
435 2225376 : return fL * std::pow((1 - psi), 1.0 / 3.0) * (1 - std::pow(psi, lambda)) +
436 2225376 : fT * std::pow(psi, 1.0 / 3.0);
437 : }
438 : }
439 :
440 : Real
441 16223040 : TriSubChannel1PhaseProblem::computeBeta(unsigned int i_gap, unsigned int iz, bool enthalpy)
442 : {
443 : auto beta = std::numeric_limits<double>::quiet_NaN();
444 16223040 : const Real & pitch = _subchannel_mesh.getPitch();
445 16223040 : const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
446 16223040 : const Real & wire_lead_length = _tri_sch_mesh.getWireLeadLength();
447 16223040 : const Real & wire_diameter = _tri_sch_mesh.getWireDiameter();
448 16223040 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
449 16223040 : auto Nr = _tri_sch_mesh._n_rings;
450 : unsigned int i_ch = chans.first;
451 : unsigned int j_ch = chans.second;
452 16223040 : auto subch_type_i = _subchannel_mesh.getSubchannelType(i_ch);
453 16223040 : auto subch_type_j = _subchannel_mesh.getSubchannelType(j_ch);
454 16223040 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
455 16223040 : auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
456 16223040 : auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
457 16223040 : auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
458 16223040 : auto Si_in = (*_S_flow_soln)(node_in_i);
459 16223040 : auto Sj_in = (*_S_flow_soln)(node_in_j);
460 16223040 : auto Si_out = (*_S_flow_soln)(node_out_i);
461 16223040 : auto Sj_out = (*_S_flow_soln)(node_out_j);
462 16223040 : auto S_total = Si_in + Sj_in + Si_out + Sj_out;
463 16223040 : auto Si = 0.5 * (Si_in + Si_out);
464 16223040 : auto Sj = 0.5 * (Sj_in + Sj_out);
465 16223040 : auto w_perim_i = 0.5 * ((*_w_perim_soln)(node_in_i) + (*_w_perim_soln)(node_out_i));
466 16223040 : auto w_perim_j = 0.5 * ((*_w_perim_soln)(node_in_j) + (*_w_perim_soln)(node_out_j));
467 16223040 : auto avg_mu = (1 / S_total) * ((*_mu_soln)(node_out_i)*Si_out + (*_mu_soln)(node_in_i)*Si_in +
468 16223040 : (*_mu_soln)(node_out_j)*Sj_out + (*_mu_soln)(node_in_j)*Sj_in);
469 16223040 : auto avg_hD = 4.0 * (Si + Sj) / (w_perim_i + w_perim_j);
470 : auto avg_massflux =
471 16223040 : 0.5 * (((*_mdot_soln)(node_in_i) + (*_mdot_soln)(node_in_j)) / (Si_in + Sj_in) +
472 16223040 : ((*_mdot_soln)(node_out_i) + (*_mdot_soln)(node_out_j)) / (Si_out + Sj_out));
473 16223040 : auto Re = avg_massflux * avg_hD / avg_mu;
474 : // crossflow area between channels i,j (dz*gap_width)
475 16223040 : auto gap = _subchannel_mesh.getGapWidth(iz, i_gap);
476 : // Calculation of flow regime
477 16223040 : auto ReL = 320.0 * std::pow(10.0, pitch / pin_diameter - 1);
478 16223040 : auto ReT = 10000.0 * std::pow(10.0, 0.7 * (pitch / pin_diameter - 1));
479 : // Calculation of Turbulent Crossflow for wire-wrapped triangular assemblies. Cheng &
480 : // Todreas (1986).
481 : // INNER SUBCHANNELS
482 16223040 : if ((subch_type_i == EChannelType::CENTER || subch_type_j == EChannelType::CENTER) &&
483 9755100 : (wire_lead_length != 0) && (wire_diameter != 0))
484 : {
485 : // Calculation of geometric parameters
486 : // wire angle
487 9165420 : auto theta = std::acos(wire_lead_length /
488 9165420 : std::sqrt(std::pow(wire_lead_length, 2) +
489 9165420 : std::pow(libMesh::pi * (pin_diameter + wire_diameter), 2)));
490 : // projected area of wire on subchannel
491 9165420 : auto Ar1 = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
492 : // bare subchannel flow area
493 : auto A1prime =
494 9165420 : (std::sqrt(3.0) / 4.0) * std::pow(pitch, 2) - libMesh::pi * std::pow(pin_diameter, 2) / 8.0;
495 : // wire-wrapped subchannel flow area
496 9165420 : auto A1 = A1prime - libMesh::pi * std::pow(wire_diameter, 2) / 8.0 / std::cos(theta);
497 : // empirical constant for mixing parameter
498 : auto Cm = 0.0;
499 : auto CmL_constant = 0.0;
500 : auto CmT_constant = 0.0;
501 :
502 9165420 : if (Nr == 1)
503 : {
504 : CmT_constant = 0.1;
505 : CmL_constant = 0.055;
506 : }
507 : else
508 : {
509 : CmT_constant = 0.14;
510 : CmL_constant = 0.077;
511 : }
512 :
513 9165420 : auto CmT = CmT_constant * std::pow((pitch - pin_diameter) / pin_diameter, -0.5);
514 9165420 : auto CmL = CmL_constant * std::pow((pitch - pin_diameter) / pin_diameter, -0.5);
515 :
516 9165420 : if (Re < ReL)
517 : {
518 : Cm = CmL;
519 : }
520 9165420 : else if (Re > ReT)
521 : {
522 : Cm = CmT;
523 : }
524 : else
525 : {
526 4261362 : auto psi = (std::log(Re) - std::log(ReL)) / (std::log(ReT) - std::log(ReL));
527 : auto gamma = 2.0 / 3.0;
528 4261362 : Cm = CmL + (CmT - CmL) * std::pow(psi, gamma);
529 : }
530 : // mixing parameter
531 9165420 : beta = Cm * std::pow(Ar1 / A1, 0.5) * std::tan(theta);
532 9165420 : }
533 : // EDGE OR CORNER SUBCHANNELS/ SWEEP FLOW
534 7057620 : else if ((subch_type_i == EChannelType::CORNER || subch_type_i == EChannelType::EDGE) &&
535 6467940 : (subch_type_j == EChannelType::CORNER || subch_type_j == EChannelType::EDGE) &&
536 6467940 : (wire_lead_length != 0) && (wire_diameter != 0))
537 : {
538 6251220 : auto theta = std::acos(wire_lead_length /
539 6251220 : std::sqrt(std::pow(wire_lead_length, 2) +
540 6251220 : std::pow(libMesh::pi * (pin_diameter + wire_diameter), 2)));
541 : // Calculation of geometric parameters
542 : // distance from pin surface to duct
543 6251220 : auto dpgap = _tri_sch_mesh.getDuctToPinGap();
544 : // Edge pitch parameter defined as pin diameter plus distance to duct wall
545 6251220 : auto w = pin_diameter + dpgap;
546 6251220 : auto Ar2 = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 4.0;
547 6251220 : auto A2prime = pitch * (w - pin_diameter / 2.0) - libMesh::pi * std::pow(pin_diameter, 2) / 8.0;
548 6251220 : auto A2 = A2prime - libMesh::pi * std::pow(wire_diameter, 2) / 8.0 / std::cos(theta);
549 : // empirical constant for mixing parameter
550 : auto Cs = 0.0;
551 : auto CsL_constant = 0.0;
552 : auto CsT_constant = 0.0;
553 6251220 : if (Nr == 1)
554 : {
555 : CsT_constant = 0.6;
556 : CsL_constant = 0.33;
557 : }
558 : else
559 : {
560 : CsT_constant = 0.75;
561 : CsL_constant = 0.413;
562 : }
563 6251220 : auto CsL = CsL_constant * std::pow(wire_lead_length / pin_diameter, 0.3);
564 6251220 : auto CsT = CsT_constant * std::pow(wire_lead_length / pin_diameter, 0.3);
565 :
566 6251220 : if (Re < ReL)
567 : {
568 : Cs = CsL;
569 : }
570 6251220 : else if (Re > ReT)
571 : {
572 : Cs = CsT;
573 : }
574 : else
575 : {
576 3178656 : auto psi = (std::log(Re) - std::log(ReL)) / (std::log(ReT) - std::log(ReL));
577 : auto gamma = 2.0 / 3.0;
578 3178656 : Cs = CsL + (CsT - CsL) * std::pow(psi, gamma);
579 : }
580 : // Calculation of turbulent mixing parameter used for sweep flow only
581 6251220 : if (enthalpy)
582 2788920 : beta = Cs * std::pow(Ar2 / A2, 0.5) * std::tan(theta);
583 : else
584 : beta = 0.0;
585 : }
586 : // Calculation of Turbulent Crossflow for bare assemblies, from Kim and Chung (2001).
587 806400 : else if ((wire_lead_length == 0) && (wire_diameter == 0))
588 : {
589 : Real gamma = 20.0; // empirical constant
590 : Real sf = 2.0 / 3.0; // shape factor
591 : Real a = 0.18;
592 : Real b = 0.2;
593 806400 : auto f = a * std::pow(Re, -b); // Rehme 1992 circular tube friction factor
594 : auto k = (1 / S_total) *
595 806400 : (_fp->k_from_p_T((*_P_soln)(node_out_i) + _P_out, (*_T_soln)(node_out_i)) * Si_out +
596 806400 : _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i)) * Si_in +
597 806400 : _fp->k_from_p_T((*_P_soln)(node_out_j) + _P_out, (*_T_soln)(node_out_j)) * Sj_out +
598 806400 : _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j)) * Sj_in);
599 : auto cp = (1 / S_total) *
600 806400 : (_fp->cp_from_p_T((*_P_soln)(node_out_i) + _P_out, (*_T_soln)(node_out_i)) * Si_out +
601 806400 : _fp->cp_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i)) * Si_in +
602 806400 : _fp->cp_from_p_T((*_P_soln)(node_out_j) + _P_out, (*_T_soln)(node_out_j)) * Sj_out +
603 806400 : _fp->cp_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j)) * Sj_in);
604 806400 : auto Pr = avg_mu * cp / k; // Prandtl number
605 806400 : auto Pr_t = Pr * (Re / gamma) * std::sqrt(f / 8.0); // Turbulent Prandtl number
606 806400 : auto delta = pitch / std::sqrt(3.0); // centroid to centroid distance
607 806400 : auto L_x = sf * delta; // axial length scale (gap is the lateral length scale)
608 806400 : auto lamda = gap / L_x; // aspect ratio
609 806400 : auto a_x = 1.0 - 2.0 * lamda * lamda / libMesh::pi; // velocity coefficient
610 806400 : auto z_FP_over_D = (2.0 * L_x / pin_diameter) *
611 806400 : (1 + (-0.5 * std::log(lamda) + 0.5 * std::log(4.0) - 0.25) * lamda * lamda);
612 806400 : auto Str = 1.0 / (0.822 * (gap / pin_diameter) + 0.144); // Strouhal number (Wu & Trupp 1994)
613 806400 : auto freq_factor = 2.0 / std::pow(gamma, 2) * std::sqrt(a / 8.0) * (avg_hD / gap);
614 806400 : auto rod_mixing = (1 / Pr_t) * lamda;
615 806400 : auto axial_mixing = a_x * z_FP_over_D * Str;
616 : // Mixing Stanton number: Stg (eq 25,Kim and Chung (2001), eq 19 (Jeong et. al 2005)
617 806400 : beta = freq_factor * (rod_mixing + axial_mixing) * std::pow(Re, -b / 2.0);
618 : }
619 : mooseAssert(beta >= 0, "beta should be positive or zero.");
620 16223040 : return beta;
621 : }
622 :
623 : Real
624 3882204 : TriSubChannel1PhaseProblem::computeAddedHeatPin(unsigned int i_ch, unsigned int iz)
625 : {
626 : // Compute axial location of nodes.
627 3882204 : auto z2 = _z_grid[iz];
628 3882204 : auto z1 = _z_grid[iz - 1];
629 3882204 : auto heated_length = _subchannel_mesh.getHeatedLength();
630 3882204 : auto unheated_length_entry = _subchannel_mesh.getHeatedLengthEntry();
631 3882204 : if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
632 2830860 : MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
633 : {
634 : // Compute the height of this element.
635 2241684 : auto dz = z2 - z1;
636 2241684 : if (_pin_mesh_exist)
637 : {
638 : double factor;
639 491904 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
640 491904 : switch (subch_type)
641 : {
642 : case EChannelType::CENTER:
643 : factor = 1.0 / 6.0;
644 : break;
645 107136 : case EChannelType::EDGE:
646 : factor = 1.0 / 4.0;
647 107136 : break;
648 : case EChannelType::CORNER:
649 : factor = 1.0 / 6.0;
650 : break;
651 : default:
652 : return 0.0; // handle invalid subch_type if needed
653 : }
654 : double heat_rate_in = 0.0;
655 : double heat_rate_out = 0.0;
656 1786752 : for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
657 : {
658 1294848 : auto * node_in = _subchannel_mesh.getPinNode(i_pin, iz - 1);
659 1294848 : auto * node_out = _subchannel_mesh.getPinNode(i_pin, iz);
660 1294848 : heat_rate_out += factor * (*_q_prime_soln)(node_out);
661 1294848 : heat_rate_in += factor * (*_q_prime_soln)(node_in);
662 : }
663 491904 : return (heat_rate_in + heat_rate_out) * dz / 2.0;
664 : }
665 : else
666 : {
667 1749780 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
668 1749780 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
669 1749780 : return ((*_q_prime_soln)(node_out) + (*_q_prime_soln)(node_in)) * dz / 2.0;
670 : }
671 : }
672 : else
673 : return 0.0;
674 : }
675 :
676 : Real
677 30240 : TriSubChannel1PhaseProblem::getSubChannelPeripheralDuctWidth(unsigned int i_ch)
678 : {
679 30240 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
680 30240 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
681 : {
682 30240 : auto width = _subchannel_mesh.getPitch();
683 30240 : if (subch_type == EChannelType::CORNER)
684 10080 : width = 2.0 / std::sqrt(3.0) *
685 10080 : (_subchannel_mesh.getPinDiameter() / 2.0 + _tri_sch_mesh.getDuctToPinGap());
686 30240 : return width;
687 : }
688 : else
689 0 : mooseError("Channel is not a perimetric subchannel ");
690 : }
691 :
692 : void
693 2262 : TriSubChannel1PhaseProblem::computeh(int iblock)
694 : {
695 2262 : unsigned int last_node = (iblock + 1) * _block_size;
696 2262 : unsigned int first_node = iblock * _block_size + 1;
697 2262 : const Real & wire_lead_length = _tri_sch_mesh.getWireLeadLength();
698 2262 : const Real & wire_diameter = _tri_sch_mesh.getWireDiameter();
699 2262 : const Real & pitch = _subchannel_mesh.getPitch();
700 2262 : const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
701 :
702 2262 : if (iblock == 0)
703 : {
704 112242 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
705 : {
706 109980 : auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
707 109980 : auto h_out = _fp->h_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
708 109980 : if (h_out < 0)
709 : {
710 0 : mooseError(
711 0 : name(), " : Calculation of negative Enthalpy h_out = : ", h_out, " Axial Level= : ", 0);
712 : }
713 109980 : _h_soln->set(node, h_out);
714 : }
715 : }
716 :
717 2262 : if (!_implicit_bool)
718 : {
719 108 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
720 : {
721 90 : auto z_grid = _subchannel_mesh.getZGrid();
722 90 : auto dz = z_grid[iz] - z_grid[iz - 1];
723 : // Calculation of average mass flux of all periphery subchannels
724 : Real edge_flux_ave = 0.0;
725 : Real mdot_sum = 0.0;
726 : Real si_sum = 0.0;
727 3870 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
728 : {
729 3780 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
730 3780 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
731 : {
732 1620 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
733 1620 : auto Si = (*_S_flow_soln)(node_in);
734 1620 : auto mdot_in = (*_mdot_soln)(node_in);
735 1620 : mdot_sum = mdot_sum + mdot_in;
736 1620 : si_sum = si_sum + Si;
737 : }
738 : }
739 90 : edge_flux_ave = mdot_sum / si_sum;
740 :
741 3870 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
742 : {
743 3780 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
744 3780 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
745 3780 : auto mdot_in = (*_mdot_soln)(node_in);
746 3780 : auto h_in = (*_h_soln)(node_in); // J/kg
747 3780 : auto volume = dz * (*_S_flow_soln)(node_in);
748 3780 : auto mdot_out = (*_mdot_soln)(node_out);
749 3780 : auto h_out = 0.0;
750 : Real sumWijh = 0.0;
751 : Real sumWijPrimeDhij = 0.0;
752 : Real sweep_enthalpy = 0.0;
753 : Real e_cond = 0.0;
754 :
755 : // Calculate added enthalpy from heatflux (Pin, Duct)
756 3780 : Real added_enthalpy = computeAddedHeatPin(i_ch, iz);
757 3780 : added_enthalpy += computeAddedHeatDuct(i_ch, iz);
758 :
759 : // Calculate net sum of enthalpy into/out-of channel i from channels j around i
760 : // (Turbulent diffusion, Diversion Crossflow, Sweep flow Enthalpy, Radial heat conduction)
761 : unsigned int counter = 0;
762 14580 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
763 : {
764 10800 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
765 10800 : auto gap = _subchannel_mesh.getGapWidth(iz, i_gap);
766 10800 : auto Sij = dz * gap;
767 : unsigned int ii_ch = chans.first; // the first subchannel next to gap i_gap
768 : unsigned int jj_ch = chans.second; // the second subchannel next to gap i_gap
769 10800 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
770 10800 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
771 10800 : auto subch_type_i = _subchannel_mesh.getSubchannelType(ii_ch);
772 10800 : auto subch_type_j = _subchannel_mesh.getSubchannelType(jj_ch);
773 : // Define donor enthalpy
774 : auto h_star = 0.0;
775 10800 : if (_Wij(i_gap, iz) > 0.0)
776 8640 : h_star = (*_h_soln)(node_in_i);
777 2160 : else if (_Wij(i_gap, iz) < 0.0)
778 2160 : h_star = (*_h_soln)(node_in_j);
779 : // Diversion crossflow
780 : // take care of the sign by applying the map, use donor cell
781 10800 : sumWijh += _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * h_star;
782 10800 : counter++;
783 : // SWEEP FLOW is calculated if i_gap is located in the periphery
784 : // and we have a wire-wrap (if i_gap is in the periphery then i_ch is in the periphery)
785 : // There are two gaps per periphery subchannel that this is true.
786 10800 : if ((subch_type_i == EChannelType::CORNER || subch_type_i == EChannelType::EDGE) &&
787 3240 : (subch_type_j == EChannelType::CORNER || subch_type_j == EChannelType::EDGE) &&
788 3240 : (wire_lead_length != 0) && (wire_diameter != 0))
789 : {
790 : // donor subchannel and node of sweep flow. The donor subchannel is the subchannel next
791 : // to i_ch that sweep flow, flows from and into i_ch
792 3240 : auto sweep_donor = _tri_sch_mesh.getSweepFlowChans(i_ch).first;
793 3240 : auto * node_sweep_donor = _subchannel_mesh.getChannelNode(sweep_donor, iz - 1);
794 : // if one of the neighbor subchannels of the periphery gap is the donor subchannel
795 : //(the other would be the i_ch) sweep enthalpy flows into i_ch
796 3240 : if ((ii_ch == sweep_donor) || (jj_ch == sweep_donor))
797 : {
798 1620 : sweep_enthalpy +=
799 1620 : computeBeta(i_gap, iz, true) * edge_flux_ave * Sij * (*_h_soln)(node_sweep_donor);
800 : }
801 : // else sweep enthalpy flows out of i_ch
802 : else
803 : {
804 1620 : sweep_enthalpy -=
805 1620 : computeBeta(i_gap, iz, true) * edge_flux_ave * Sij * (*_h_soln)(node_in);
806 : }
807 : }
808 : // Inner gap
809 : // Turbulent Diffusion
810 : else
811 : {
812 7560 : sumWijPrimeDhij +=
813 7560 : _WijPrime(i_gap, iz) * (2 * h_in - (*_h_soln)(node_in_j) - (*_h_soln)(node_in_i));
814 : }
815 :
816 : // compute the radial heat conduction through the gaps
817 : Real dist_ij = pitch;
818 :
819 10800 : if (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::EDGE)
820 : {
821 1080 : dist_ij = pitch;
822 : }
823 9720 : else if ((subch_type_i == EChannelType::CORNER && subch_type_j == EChannelType::EDGE) ||
824 9540 : (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::CORNER))
825 : {
826 2160 : dist_ij = pitch;
827 : }
828 : else
829 : {
830 7560 : dist_ij = pitch / std::sqrt(3);
831 : }
832 :
833 10800 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
834 10800 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
835 : auto shape_factor =
836 10800 : 0.66 * (pitch / pin_diameter) *
837 10800 : std::pow((_subchannel_mesh.getGapWidth(iz, i_gap) / pin_diameter), -0.3);
838 10800 : if (ii_ch == i_ch)
839 : {
840 5400 : e_cond += 0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
841 5400 : ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
842 : }
843 : else
844 : {
845 5400 : e_cond += -0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
846 5400 : ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
847 : }
848 : }
849 :
850 : // compute the axial heat conduction between current and lower axial node
851 3780 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz);
852 3780 : auto * node_in_j = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
853 3780 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
854 3780 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
855 3780 : auto Si = (*_S_flow_soln)(node_in_i);
856 3780 : auto dist_ij = z_grid[iz] - z_grid[iz - 1];
857 :
858 3780 : e_cond += 0.5 * (thcon_i + thcon_j) * Si * ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) /
859 : dist_ij;
860 :
861 3780 : unsigned int nz = _subchannel_mesh.getNumOfAxialCells();
862 : // compute the axial heat conduction between current and upper axial node
863 3780 : if (iz < nz)
864 : {
865 3024 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz);
866 3024 : auto * node_in_j = _subchannel_mesh.getChannelNode(i_ch, iz + 1);
867 3024 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
868 3024 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
869 3024 : auto Si = (*_S_flow_soln)(node_in_i);
870 3024 : auto dist_ij = z_grid[iz + 1] - z_grid[iz];
871 3024 : e_cond += 0.5 * (thcon_i + thcon_j) * Si *
872 3024 : ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
873 : }
874 :
875 : // end of radial heat conduction calc.
876 3780 : h_out =
877 7560 : (mdot_in * h_in - sumWijh - sumWijPrimeDhij + added_enthalpy + e_cond + sweep_enthalpy +
878 3780 : _TR * _rho_soln->old(node_out) * _h_soln->old(node_out) * volume / _dt) /
879 3780 : (mdot_out + _TR * (*_rho_soln)(node_out)*volume / _dt);
880 3780 : if (h_out < 0)
881 : {
882 0 : mooseError(name(),
883 : " : Calculation of negative Enthalpy h_out = : ",
884 : h_out,
885 : " Axial Level= : ",
886 : iz);
887 : }
888 3780 : _h_soln->set(node_out, h_out); // J/kg
889 : }
890 90 : }
891 : }
892 : else
893 : {
894 2244 : LibmeshPetscCall(MatZeroEntries(_hc_time_derivative_mat));
895 2244 : LibmeshPetscCall(MatZeroEntries(_hc_advective_derivative_mat));
896 2244 : LibmeshPetscCall(MatZeroEntries(_hc_cross_derivative_mat));
897 2244 : LibmeshPetscCall(MatZeroEntries(_hc_axial_heat_conduction_mat));
898 2244 : LibmeshPetscCall(MatZeroEntries(_hc_radial_heat_conduction_mat));
899 2244 : LibmeshPetscCall(MatZeroEntries(_hc_sweep_enthalpy_mat));
900 :
901 2244 : LibmeshPetscCall(VecZeroEntries(_hc_time_derivative_rhs));
902 2244 : LibmeshPetscCall(VecZeroEntries(_hc_advective_derivative_rhs));
903 2244 : LibmeshPetscCall(VecZeroEntries(_hc_cross_derivative_rhs));
904 2244 : LibmeshPetscCall(VecZeroEntries(_hc_added_heat_rhs));
905 2244 : LibmeshPetscCall(VecZeroEntries(_hc_axial_heat_conduction_rhs));
906 2244 : LibmeshPetscCall(VecZeroEntries(_hc_radial_heat_conduction_rhs));
907 2244 : LibmeshPetscCall(VecZeroEntries(_hc_sweep_enthalpy_rhs));
908 :
909 2244 : LibmeshPetscCall(MatZeroEntries(_hc_sys_h_mat));
910 2244 : LibmeshPetscCall(VecZeroEntries(_hc_sys_h_rhs));
911 :
912 80844 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
913 : {
914 78600 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
915 78600 : auto pitch = _subchannel_mesh.getPitch();
916 78600 : auto pin_diameter = _subchannel_mesh.getPinDiameter();
917 78600 : auto iz_ind = iz - first_node;
918 :
919 3921240 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
920 : {
921 3842640 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
922 3842640 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
923 3842640 : auto S_in = (*_S_flow_soln)(node_in);
924 3842640 : auto S_out = (*_S_flow_soln)(node_out);
925 3842640 : auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
926 3842640 : auto volume = dz * S_interp;
927 :
928 : PetscScalar Pe = 0.5;
929 3842640 : if (_interpolation_scheme == 3)
930 : {
931 : // Compute the Peclet number
932 942984 : auto w_perim_in = (*_w_perim_soln)(node_in);
933 942984 : auto w_perim_out = (*_w_perim_soln)(node_out);
934 942984 : auto w_perim_interp = this->computeInterpolatedValue(w_perim_out, w_perim_in, 0.5);
935 942984 : auto K_in = _fp->k_from_p_T((*_P_soln)(node_in) + _P_out, (*_T_soln)(node_in));
936 942984 : auto K_out = _fp->k_from_p_T((*_P_soln)(node_out) + _P_out, (*_T_soln)(node_out));
937 942984 : auto K = this->computeInterpolatedValue(K_out, K_in, 0.5);
938 942984 : auto cp_in = _fp->cp_from_p_T((*_P_soln)(node_in) + _P_out, (*_T_soln)(node_in));
939 942984 : auto cp_out = _fp->cp_from_p_T((*_P_soln)(node_out) + _P_out, (*_T_soln)(node_out));
940 942984 : auto cp = this->computeInterpolatedValue(cp_out, cp_in, 0.5);
941 : auto mdot_loc =
942 942984 : this->computeInterpolatedValue((*_mdot_soln)(node_out), (*_mdot_soln)(node_in), 0.5);
943 : // hydraulic diameter in the i direction
944 942984 : auto Dh_i = 4.0 * S_interp / w_perim_interp;
945 942984 : Pe = mdot_loc * Dh_i * cp / (K * S_interp) * (mdot_loc / std::abs(mdot_loc));
946 : }
947 3842640 : auto alpha = computeInterpolationCoefficients(Pe);
948 :
949 : // Time derivative term
950 3842640 : if (iz == first_node)
951 : {
952 : PetscScalar value_vec_tt =
953 109224 : -1.0 * _TR * alpha * (*_rho_soln)(node_in) * (*_h_soln)(node_in)*volume / _dt;
954 109224 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
955 109224 : LibmeshPetscCall(
956 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
957 : }
958 : else
959 : {
960 3733416 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
961 3733416 : PetscInt col_tt = i_ch + _n_channels * (iz_ind - 1);
962 3733416 : PetscScalar value_tt = _TR * alpha * (*_rho_soln)(node_in)*volume / _dt;
963 3733416 : LibmeshPetscCall(MatSetValues(
964 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
965 : }
966 :
967 : // Adding diagonal elements
968 3842640 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
969 3842640 : PetscInt col_tt = i_ch + _n_channels * iz_ind;
970 3842640 : PetscScalar value_tt = _TR * (1.0 - alpha) * (*_rho_soln)(node_out)*volume / _dt;
971 3842640 : LibmeshPetscCall(MatSetValues(
972 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
973 :
974 : // Adding RHS elements
975 : PetscScalar rho_old_interp =
976 3842640 : computeInterpolatedValue(_rho_soln->old(node_out), _rho_soln->old(node_in), Pe);
977 : PetscScalar h_old_interp =
978 3842640 : computeInterpolatedValue(_h_soln->old(node_out), _h_soln->old(node_in), Pe);
979 3842640 : PetscScalar value_vec_tt = _TR * rho_old_interp * h_old_interp * volume / _dt;
980 3842640 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
981 3842640 : LibmeshPetscCall(
982 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
983 :
984 : // Advective derivative term
985 3842640 : if (iz == first_node)
986 : {
987 109224 : PetscInt row_at = i_ch + _n_channels * iz_ind;
988 109224 : PetscScalar value_at = alpha * (*_mdot_soln)(node_in) * (*_h_soln)(node_in);
989 109224 : LibmeshPetscCall(
990 : VecSetValues(_hc_advective_derivative_rhs, 1, &row_at, &value_at, ADD_VALUES));
991 :
992 109224 : value_at = alpha * (*_mdot_soln)(node_out) - (1 - alpha) * (*_mdot_soln)(node_in);
993 109224 : PetscInt col_at = i_ch + _n_channels * iz_ind;
994 109224 : LibmeshPetscCall(MatSetValues(
995 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
996 :
997 109224 : value_at = (1 - alpha) * (*_mdot_soln)(node_out);
998 109224 : col_at = i_ch + _n_channels * (iz_ind + 1);
999 109224 : LibmeshPetscCall(MatSetValues(
1000 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
1001 : }
1002 3733416 : else if (iz == last_node)
1003 : {
1004 109224 : PetscInt row_at = i_ch + _n_channels * iz_ind;
1005 109224 : PetscScalar value_at = 1.0 * (*_mdot_soln)(node_out);
1006 109224 : PetscInt col_at = i_ch + _n_channels * iz_ind;
1007 109224 : LibmeshPetscCall(MatSetValues(
1008 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
1009 :
1010 109224 : value_at = -1.0 * (*_mdot_soln)(node_in);
1011 109224 : col_at = i_ch + _n_channels * (iz_ind - 1);
1012 109224 : LibmeshPetscCall(MatSetValues(
1013 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
1014 : }
1015 : else
1016 : {
1017 3624192 : PetscInt row_at = i_ch + _n_channels * iz_ind;
1018 : PetscInt col_at;
1019 :
1020 3624192 : PetscScalar value_at = -alpha * (*_mdot_soln)(node_in);
1021 3624192 : col_at = i_ch + _n_channels * (iz_ind - 1);
1022 3624192 : LibmeshPetscCall(MatSetValues(
1023 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
1024 :
1025 3624192 : value_at = alpha * (*_mdot_soln)(node_out) - (1 - alpha) * (*_mdot_soln)(node_in);
1026 3624192 : col_at = i_ch + _n_channels * iz_ind;
1027 3624192 : LibmeshPetscCall(MatSetValues(
1028 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
1029 :
1030 3624192 : value_at = (1 - alpha) * (*_mdot_soln)(node_out);
1031 3624192 : col_at = i_ch + _n_channels * (iz_ind + 1);
1032 3624192 : LibmeshPetscCall(MatSetValues(
1033 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
1034 : }
1035 :
1036 : // Axial heat conduction
1037 3842640 : auto * node_center = _subchannel_mesh.getChannelNode(i_ch, iz);
1038 3842640 : auto K_center = _fp->k_from_p_T((*_P_soln)(node_center) + _P_out, (*_T_soln)(node_center));
1039 : auto cp_center =
1040 3842640 : _fp->cp_from_p_T((*_P_soln)(node_center) + _P_out, (*_T_soln)(node_center));
1041 3842640 : auto diff_center = K_center / (cp_center + 1e-15);
1042 :
1043 3842640 : if (iz == first_node)
1044 : {
1045 109224 : auto * node_top = _subchannel_mesh.getChannelNode(i_ch, iz + 1);
1046 109224 : auto * node_bottom = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1047 : auto K_bottom =
1048 109224 : _fp->k_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1049 109224 : auto K_top = _fp->k_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
1050 : auto cp_bottom =
1051 109224 : _fp->cp_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1052 109224 : auto cp_top = _fp->cp_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
1053 109224 : auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
1054 109224 : auto diff_top = K_top / (cp_top + 1e-15);
1055 :
1056 109224 : auto dz_up = _z_grid[iz + 1] - _z_grid[iz];
1057 109224 : auto dz_down = _z_grid[iz] - _z_grid[iz - 1];
1058 : auto S_up =
1059 109224 : computeInterpolatedValue((*_S_flow_soln)(node_top), (*_S_flow_soln)(node_center));
1060 : auto S_down =
1061 109224 : computeInterpolatedValue((*_S_flow_soln)(node_center), (*_S_flow_soln)(node_bottom));
1062 109224 : auto diff_up = computeInterpolatedValue(diff_top, diff_center);
1063 109224 : auto diff_down = computeInterpolatedValue(diff_center, diff_bottom);
1064 :
1065 : // Diagonal value
1066 109224 : PetscInt row_at = i_ch + _n_channels * iz_ind;
1067 109224 : PetscInt col_at = i_ch + _n_channels * iz_ind;
1068 109224 : PetscScalar value_at = diff_up * S_up / dz_up + diff_down * S_down / dz_down;
1069 109224 : LibmeshPetscCall(MatSetValues(
1070 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1071 :
1072 : // Bottom value
1073 109224 : value_at = 1.0 * diff_down * S_down / dz_down * (*_h_soln)(node_bottom);
1074 109224 : LibmeshPetscCall(
1075 : VecSetValues(_hc_axial_heat_conduction_rhs, 1, &row_at, &value_at, ADD_VALUES));
1076 :
1077 : // Top value
1078 109224 : col_at = i_ch + _n_channels * (iz_ind + 1);
1079 109224 : value_at = -diff_up * S_up / dz_up;
1080 109224 : LibmeshPetscCall(MatSetValues(
1081 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1082 : }
1083 3733416 : else if (iz == last_node)
1084 : {
1085 109224 : auto * node_bottom = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1086 : auto K_bottom =
1087 109224 : _fp->k_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1088 : auto cp_bottom =
1089 109224 : _fp->cp_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1090 109224 : auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
1091 :
1092 109224 : auto dz_down = _z_grid[iz] - _z_grid[iz - 1];
1093 109224 : auto S_down = 0.5 * ((*_S_flow_soln)(node_center) + (*_S_flow_soln)(node_bottom));
1094 109224 : auto diff_down = 0.5 * (diff_center + diff_bottom);
1095 :
1096 : // Diagonal value
1097 109224 : PetscInt row_at = i_ch + _n_channels * iz_ind;
1098 109224 : PetscInt col_at = i_ch + _n_channels * iz_ind;
1099 109224 : PetscScalar value_at = diff_down * S_down / dz_down;
1100 109224 : LibmeshPetscCall(MatSetValues(
1101 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1102 :
1103 : // Bottom value
1104 109224 : col_at = i_ch + _n_channels * (iz_ind - 1);
1105 109224 : value_at = -diff_down * S_down / dz_down;
1106 109224 : LibmeshPetscCall(MatSetValues(
1107 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1108 :
1109 : // Outflow derivative
1110 : /// TODO: Current axial derivative is zero - check if outflow conditions may make a difference
1111 : // value_at = -1.0 * (*_mdot_soln)(node_center) * (*_h_soln)(node_center);
1112 : // VecSetValues(_hc_axial_heat_conduction_rhs, 1, &row_at, &value_at, ADD_VALUES);
1113 : }
1114 : else
1115 : {
1116 3624192 : auto * node_top = _subchannel_mesh.getChannelNode(i_ch, iz + 1);
1117 3624192 : auto * node_bottom = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1118 : auto K_bottom =
1119 3624192 : _fp->k_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1120 3624192 : auto K_top = _fp->k_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
1121 : auto cp_bottom =
1122 3624192 : _fp->cp_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1123 3624192 : auto cp_top = _fp->cp_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
1124 3624192 : auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
1125 3624192 : auto diff_top = K_top / (cp_top + 1e-15);
1126 :
1127 3624192 : auto dz_up = _z_grid[iz + 1] - _z_grid[iz];
1128 3624192 : auto dz_down = _z_grid[iz] - _z_grid[iz - 1];
1129 : auto S_up =
1130 3624192 : computeInterpolatedValue((*_S_flow_soln)(node_top), (*_S_flow_soln)(node_center));
1131 : auto S_down =
1132 3624192 : computeInterpolatedValue((*_S_flow_soln)(node_center), (*_S_flow_soln)(node_bottom));
1133 3624192 : auto diff_up = computeInterpolatedValue(diff_top, diff_center);
1134 3624192 : auto diff_down = computeInterpolatedValue(diff_center, diff_bottom);
1135 :
1136 : // Diagonal value
1137 3624192 : PetscInt row_at = i_ch + _n_channels * iz_ind;
1138 3624192 : PetscInt col_at = i_ch + _n_channels * iz_ind;
1139 3624192 : PetscScalar value_at = diff_up * S_up / dz_up + diff_down * S_down / dz_down;
1140 3624192 : LibmeshPetscCall(MatSetValues(
1141 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1142 :
1143 : // Bottom value
1144 3624192 : col_at = i_ch + _n_channels * (iz_ind - 1);
1145 3624192 : value_at = -diff_down * S_down / dz_down;
1146 3624192 : LibmeshPetscCall(MatSetValues(
1147 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1148 :
1149 : // Top value
1150 3624192 : col_at = i_ch + _n_channels * (iz_ind + 1);
1151 3624192 : value_at = -diff_up * S_up / dz_up;
1152 3624192 : LibmeshPetscCall(MatSetValues(
1153 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1154 : }
1155 :
1156 : // Radial Terms
1157 : unsigned int counter = 0;
1158 : unsigned int cross_index = iz;
1159 : // Real radial_heat_conduction(0.0);
1160 14898960 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
1161 : {
1162 11056320 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
1163 : unsigned int ii_ch = chans.first;
1164 : unsigned int jj_ch = chans.second;
1165 11056320 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
1166 11056320 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
1167 : PetscScalar h_star;
1168 : // figure out donor axial velocity
1169 11056320 : if (_Wij(i_gap, cross_index) > 0.0)
1170 : {
1171 6766044 : if (iz == first_node)
1172 : {
1173 208632 : h_star = (*_h_soln)(node_in_i);
1174 417264 : PetscScalar value_vec_ct = -1.0 * alpha *
1175 208632 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1176 208632 : _Wij(i_gap, cross_index) * h_star;
1177 208632 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
1178 208632 : LibmeshPetscCall(VecSetValues(
1179 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
1180 : }
1181 : else
1182 : {
1183 6557412 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1184 6557412 : _Wij(i_gap, cross_index);
1185 6557412 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1186 6557412 : PetscInt col_ct = ii_ch + _n_channels * (iz_ind - 1);
1187 6557412 : LibmeshPetscCall(MatSetValues(
1188 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
1189 : }
1190 6766044 : PetscScalar value_ct = (1.0 - alpha) *
1191 6766044 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1192 6766044 : _Wij(i_gap, cross_index);
1193 6766044 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1194 6766044 : PetscInt col_ct = ii_ch + _n_channels * iz_ind;
1195 6766044 : LibmeshPetscCall(MatSetValues(
1196 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
1197 : }
1198 4290276 : else if (_Wij(i_gap, cross_index) < 0.0) // _Wij=0 operations not necessary
1199 : {
1200 3941184 : if (iz == first_node)
1201 : {
1202 103980 : h_star = (*_h_soln)(node_in_j);
1203 207960 : PetscScalar value_vec_ct = -1.0 * alpha *
1204 103980 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1205 103980 : _Wij(i_gap, cross_index) * h_star;
1206 103980 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
1207 103980 : LibmeshPetscCall(VecSetValues(
1208 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
1209 : }
1210 : else
1211 : {
1212 3837204 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1213 3837204 : _Wij(i_gap, cross_index);
1214 3837204 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1215 3837204 : PetscInt col_ct = jj_ch + _n_channels * (iz_ind - 1);
1216 3837204 : LibmeshPetscCall(MatSetValues(
1217 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
1218 : }
1219 3941184 : PetscScalar value_ct = (1.0 - alpha) *
1220 3941184 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1221 3941184 : _Wij(i_gap, cross_index);
1222 3941184 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1223 3941184 : PetscInt col_ct = jj_ch + _n_channels * iz_ind;
1224 3941184 : LibmeshPetscCall(MatSetValues(
1225 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
1226 : }
1227 :
1228 : // Turbulent cross flows
1229 11056320 : if (iz == first_node)
1230 : {
1231 : PetscScalar value_vec_ct =
1232 314208 : -2.0 * alpha * (*_h_soln)(node_in)*_WijPrime(i_gap, cross_index);
1233 314208 : value_vec_ct += alpha * (*_h_soln)(node_in_j)*_WijPrime(i_gap, cross_index);
1234 314208 : value_vec_ct += alpha * (*_h_soln)(node_in_i)*_WijPrime(i_gap, cross_index);
1235 314208 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
1236 314208 : LibmeshPetscCall(
1237 : VecSetValues(_hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
1238 : }
1239 : else
1240 : {
1241 10742112 : PetscScalar value_center_ct = 2.0 * alpha * _WijPrime(i_gap, cross_index);
1242 10742112 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1243 10742112 : PetscInt col_ct = i_ch + _n_channels * (iz_ind - 1);
1244 10742112 : LibmeshPetscCall(MatSetValues(
1245 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
1246 :
1247 10742112 : PetscScalar value_left_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
1248 10742112 : row_ct = i_ch + _n_channels * iz_ind;
1249 10742112 : col_ct = jj_ch + _n_channels * (iz_ind - 1);
1250 10742112 : LibmeshPetscCall(MatSetValues(
1251 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
1252 :
1253 10742112 : PetscScalar value_right_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
1254 10742112 : row_ct = i_ch + _n_channels * iz_ind;
1255 10742112 : col_ct = ii_ch + _n_channels * (iz_ind - 1);
1256 10742112 : LibmeshPetscCall(MatSetValues(
1257 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
1258 : }
1259 11056320 : PetscScalar value_center_ct = 2.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
1260 11056320 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1261 11056320 : PetscInt col_ct = i_ch + _n_channels * iz_ind;
1262 11056320 : LibmeshPetscCall(MatSetValues(
1263 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
1264 :
1265 11056320 : PetscScalar value_left_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
1266 11056320 : row_ct = i_ch + _n_channels * iz_ind;
1267 11056320 : col_ct = jj_ch + _n_channels * iz_ind;
1268 11056320 : LibmeshPetscCall(MatSetValues(
1269 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
1270 :
1271 11056320 : PetscScalar value_right_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
1272 11056320 : row_ct = i_ch + _n_channels * iz_ind;
1273 11056320 : col_ct = ii_ch + _n_channels * iz_ind;
1274 11056320 : LibmeshPetscCall(MatSetValues(
1275 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
1276 :
1277 : // Radial heat conduction
1278 11056320 : auto subch_type_i = _subchannel_mesh.getSubchannelType(ii_ch);
1279 11056320 : auto subch_type_j = _subchannel_mesh.getSubchannelType(jj_ch);
1280 : Real dist_ij = pitch;
1281 :
1282 11056320 : if (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::EDGE)
1283 : {
1284 : dist_ij = pitch;
1285 : }
1286 9951840 : else if ((subch_type_i == EChannelType::CORNER && subch_type_j == EChannelType::EDGE) ||
1287 9794640 : (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::CORNER))
1288 : {
1289 : dist_ij = pitch;
1290 : }
1291 : else
1292 : {
1293 8065440 : dist_ij = pitch / std::sqrt(3);
1294 : }
1295 :
1296 11056320 : auto Sij = dz * _subchannel_mesh.getGapWidth(iz, i_gap);
1297 11056320 : auto K_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
1298 11056320 : auto K_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
1299 11056320 : auto cp_i = _fp->cp_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
1300 11056320 : auto cp_j = _fp->cp_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
1301 11056320 : auto A_i = K_i / cp_i;
1302 11056320 : auto A_j = K_j / cp_j;
1303 11056320 : auto harm_A = 2.0 * A_i * A_j / (A_i + A_j);
1304 : auto shape_factor =
1305 11056320 : 0.66 * (pitch / pin_diameter) *
1306 11056320 : std::pow((_subchannel_mesh.getGapWidth(iz, i_gap) / pin_diameter), -0.3);
1307 : // auto base_value = 0.5 * (A_i + A_j) * Sij * shape_factor / dist_ij;
1308 11056320 : auto base_value = harm_A * shape_factor * Sij / dist_ij;
1309 11056320 : auto neg_base_value = -1.0 * base_value;
1310 :
1311 11056320 : row_ct = ii_ch + _n_channels * iz_ind;
1312 11056320 : col_ct = ii_ch + _n_channels * iz_ind;
1313 11056320 : LibmeshPetscCall(MatSetValues(
1314 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &base_value, ADD_VALUES));
1315 :
1316 11056320 : row_ct = jj_ch + _n_channels * iz_ind;
1317 11056320 : col_ct = jj_ch + _n_channels * iz_ind;
1318 11056320 : LibmeshPetscCall(MatSetValues(
1319 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &base_value, ADD_VALUES));
1320 :
1321 11056320 : row_ct = ii_ch + _n_channels * iz_ind;
1322 11056320 : col_ct = jj_ch + _n_channels * iz_ind;
1323 11056320 : LibmeshPetscCall(MatSetValues(
1324 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &neg_base_value, ADD_VALUES));
1325 :
1326 11056320 : row_ct = jj_ch + _n_channels * iz_ind;
1327 11056320 : col_ct = ii_ch + _n_channels * iz_ind;
1328 11056320 : LibmeshPetscCall(MatSetValues(
1329 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &neg_base_value, ADD_VALUES));
1330 11056320 : counter++;
1331 : }
1332 :
1333 : // Compute the sweep flow enthalpy change
1334 : // Calculation of average mass flux of all periphery subchannels
1335 : Real edge_flux_ave = 0.0;
1336 : Real mdot_sum = 0.0;
1337 : Real si_sum = 0.0;
1338 226819440 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1339 : {
1340 222976800 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
1341 222976800 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
1342 : {
1343 78222240 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1344 78222240 : auto Si = (*_S_flow_soln)(node_in);
1345 78222240 : auto mdot_in = (*_mdot_soln)(node_in);
1346 78222240 : mdot_sum = mdot_sum + mdot_in;
1347 78222240 : si_sum = si_sum + Si;
1348 : }
1349 : }
1350 3842640 : edge_flux_ave = mdot_sum / si_sum;
1351 3842640 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
1352 : PetscScalar sweep_enthalpy = 0.0;
1353 3842640 : if ((subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER) &&
1354 1495440 : (wire_diameter != 0.0) && (wire_lead_length != 0.0))
1355 : {
1356 : auto beta_in = std::numeric_limits<double>::quiet_NaN();
1357 : auto beta_out = std::numeric_limits<double>::quiet_NaN();
1358 : // donor sweep channel for i_ch
1359 1392840 : auto sweep_donor = _tri_sch_mesh.getSweepFlowChans(i_ch).first;
1360 1392840 : auto * node_sweep_donor = _subchannel_mesh.getChannelNode(sweep_donor, iz - 1);
1361 : // Calculation of turbulent mixing parameter
1362 5133960 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
1363 : {
1364 3741120 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
1365 : unsigned int ii_ch = chans.first;
1366 : unsigned int jj_ch = chans.second;
1367 3741120 : auto subch_type_i = _subchannel_mesh.getSubchannelType(ii_ch);
1368 3741120 : auto subch_type_j = _subchannel_mesh.getSubchannelType(jj_ch);
1369 3741120 : if ((subch_type_i == EChannelType::CORNER || subch_type_i == EChannelType::EDGE) &&
1370 2785680 : (subch_type_j == EChannelType::CORNER || subch_type_j == EChannelType::EDGE))
1371 : {
1372 2785680 : if ((ii_ch == sweep_donor) || (jj_ch == sweep_donor))
1373 : {
1374 1392840 : beta_in = computeBeta(i_gap, iz, true);
1375 : }
1376 : else
1377 : {
1378 1392840 : beta_out = computeBeta(i_gap, iz, true);
1379 : }
1380 : }
1381 : }
1382 : // Abort execution if required values are unset
1383 : mooseAssert(!std::isnan(beta_in),
1384 : "beta_in was not set. Check gap logic for i_ch = " + std::to_string(i_ch) +
1385 : ", iz = " + std::to_string(iz));
1386 : mooseAssert(!std::isnan(beta_out),
1387 : "beta_out was not set. Check gap logic for i_ch = " + std::to_string(i_ch) +
1388 : ", iz = " + std::to_string(iz));
1389 :
1390 1392840 : auto gap = _tri_sch_mesh.getDuctToPinGap();
1391 1392840 : auto Sij = dz * gap;
1392 1392840 : auto wsweep_in = edge_flux_ave * beta_in * Sij;
1393 1392840 : auto wsweep_out = edge_flux_ave * beta_out * Sij;
1394 1392840 : auto sweep_hin = (*_h_soln)(node_sweep_donor);
1395 1392840 : auto sweep_hout = (*_h_soln)(node_in);
1396 1392840 : sweep_enthalpy = (wsweep_in * sweep_hin - wsweep_out * sweep_hout);
1397 :
1398 1392840 : if (iz == first_node)
1399 : {
1400 40644 : PetscInt row_sh = i_ch + _n_channels * iz_ind;
1401 40644 : PetscScalar value_hs = -sweep_enthalpy;
1402 40644 : LibmeshPetscCall(
1403 : VecSetValues(_hc_sweep_enthalpy_rhs, 1, &row_sh, &value_hs, ADD_VALUES));
1404 : }
1405 : else
1406 : {
1407 : // coefficient of sweep_hin
1408 1352196 : PetscInt row_sh = i_ch + _n_channels * (iz_ind - 1);
1409 1352196 : PetscInt col_sh = i_ch + _n_channels * (iz_ind - 1);
1410 1352196 : LibmeshPetscCall(MatSetValues(
1411 : _hc_sweep_enthalpy_mat, 1, &row_sh, 1, &col_sh, &wsweep_out, ADD_VALUES));
1412 1352196 : PetscInt col_sh_l = sweep_donor + _n_channels * (iz_ind - 1);
1413 1352196 : PetscScalar neg_sweep_in = -1.0 * wsweep_in;
1414 : // coefficient of sweep_hout
1415 1352196 : LibmeshPetscCall(MatSetValues(
1416 : _hc_sweep_enthalpy_mat, 1, &row_sh, 1, &col_sh_l, &(neg_sweep_in), ADD_VALUES));
1417 : }
1418 : }
1419 :
1420 : // Add heat enthalpy from pin and/or duct
1421 3842640 : PetscScalar added_enthalpy = computeAddedHeatPin(i_ch, iz);
1422 3842640 : added_enthalpy += computeAddedHeatDuct(i_ch, iz);
1423 3842640 : PetscInt row_vec_ht = i_ch + _n_channels * iz_ind;
1424 3842640 : LibmeshPetscCall(
1425 : VecSetValues(_hc_added_heat_rhs, 1, &row_vec_ht, &added_enthalpy, ADD_VALUES));
1426 : }
1427 : }
1428 : // Assembling system
1429 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
1430 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
1431 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
1432 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
1433 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
1434 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
1435 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_axial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1436 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_axial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1437 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_radial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1438 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_radial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1439 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sweep_enthalpy_mat, MAT_FINAL_ASSEMBLY));
1440 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sweep_enthalpy_mat, MAT_FINAL_ASSEMBLY));
1441 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1442 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1443 : // Add all matrices together
1444 2244 : LibmeshPetscCall(
1445 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1446 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1447 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1448 2244 : LibmeshPetscCall(
1449 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1450 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1451 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1452 2244 : LibmeshPetscCall(
1453 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1454 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1455 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1456 2244 : LibmeshPetscCall(
1457 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_axial_heat_conduction_mat, DIFFERENT_NONZERO_PATTERN));
1458 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1459 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1460 2244 : LibmeshPetscCall(
1461 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_radial_heat_conduction_mat, DIFFERENT_NONZERO_PATTERN));
1462 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1463 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1464 2244 : LibmeshPetscCall(
1465 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_sweep_enthalpy_mat, DIFFERENT_NONZERO_PATTERN));
1466 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1467 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1468 2244 : if (_verbose_subchannel)
1469 1758 : _console << "Block: " << iblock << " - Enthalpy conservation matrix assembled" << std::endl;
1470 : // RHS
1471 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_time_derivative_rhs));
1472 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_advective_derivative_rhs));
1473 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_cross_derivative_rhs));
1474 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_added_heat_rhs));
1475 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_axial_heat_conduction_rhs));
1476 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_radial_heat_conduction_rhs));
1477 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_sweep_enthalpy_rhs));
1478 :
1479 2244 : if (_segregated_bool || (!_monolithic_thermal_bool))
1480 : {
1481 : // Assembly the matrix system
1482 : KSP ksploc;
1483 : PC pc;
1484 : Vec sol;
1485 1380 : LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &sol));
1486 1380 : LibmeshPetscCall(KSPCreate(PETSC_COMM_SELF, &ksploc));
1487 1380 : LibmeshPetscCall(KSPSetOperators(ksploc, _hc_sys_h_mat, _hc_sys_h_mat));
1488 1380 : LibmeshPetscCall(KSPGetPC(ksploc, &pc));
1489 1380 : LibmeshPetscCall(PCSetType(pc, PCJACOBI));
1490 1380 : LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit));
1491 1380 : LibmeshPetscCall(KSPSetOptionsPrefix(ksploc, "h_sys_"));
1492 1380 : LibmeshPetscCall(KSPSetFromOptions(ksploc));
1493 1380 : LibmeshPetscCall(KSPSolve(ksploc, _hc_sys_h_rhs, sol));
1494 : // VecView(sol, PETSC_VIEWER_STDOUT_WORLD);
1495 : PetscScalar * xx;
1496 1380 : LibmeshPetscCall(VecGetArray(sol, &xx));
1497 51300 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1498 : {
1499 49920 : auto iz_ind = iz - first_node;
1500 2688000 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1501 : {
1502 2638080 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
1503 2638080 : auto h_out = xx[iz_ind * _n_channels + i_ch];
1504 2638080 : if (h_out < 0)
1505 : {
1506 0 : mooseError(name(),
1507 : " : Calculation of negative Enthalpy h_out = : ",
1508 : h_out,
1509 : " Axial Level= : ",
1510 : iz);
1511 : }
1512 2638080 : _h_soln->set(node_out, h_out);
1513 : }
1514 : }
1515 1380 : LibmeshPetscCall(KSPDestroy(&ksploc));
1516 1380 : LibmeshPetscCall(VecDestroy(&sol));
1517 : }
1518 : }
1519 2262 : }
|
