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 204 : TriSubChannel1PhaseProblem::validParams()
21 : {
22 204 : InputParameters params = SubChannel1PhaseProblem::validParams();
23 204 : params.addClassDescription("Solver class for subchannels in a triangular lattice assembly and "
24 : "bare/wire-wrapped fuel pins");
25 204 : return params;
26 0 : }
27 :
28 102 : TriSubChannel1PhaseProblem::TriSubChannel1PhaseProblem(const InputParameters & params)
29 : : SubChannel1PhaseProblem(params),
30 102 : _tri_sch_mesh(SCM::getMesh<TriSubChannelMesh>(_subchannel_mesh))
31 : {
32 : // Initializing heat conduction system
33 102 : LibmeshPetscCall(createPetscMatrix(
34 : _hc_axial_heat_conduction_mat, _block_size * _n_channels, _block_size * _n_channels));
35 102 : LibmeshPetscCall(createPetscVector(_hc_axial_heat_conduction_rhs, _block_size * _n_channels));
36 102 : LibmeshPetscCall(createPetscMatrix(
37 : _hc_radial_heat_conduction_mat, _block_size * _n_channels, _block_size * _n_channels));
38 102 : LibmeshPetscCall(createPetscVector(_hc_radial_heat_conduction_rhs, _block_size * _n_channels));
39 102 : LibmeshPetscCall(createPetscMatrix(
40 : _hc_sweep_enthalpy_mat, _block_size * _n_channels, _block_size * _n_channels));
41 102 : LibmeshPetscCall(createPetscVector(_hc_sweep_enthalpy_rhs, _block_size * _n_channels));
42 102 : }
43 :
44 204 : TriSubChannel1PhaseProblem::~TriSubChannel1PhaseProblem()
45 : {
46 102 : PetscErrorCode ierr = cleanUp();
47 102 : if (ierr)
48 0 : mooseError(name(), ": Error in memory cleanup");
49 204 : }
50 :
51 : PetscErrorCode
52 102 : TriSubChannel1PhaseProblem::cleanUp()
53 : {
54 : PetscFunctionBegin;
55 : // Clean up heat conduction system
56 102 : LibmeshPetscCall(MatDestroy(&_hc_axial_heat_conduction_mat));
57 102 : LibmeshPetscCall(VecDestroy(&_hc_axial_heat_conduction_rhs));
58 102 : LibmeshPetscCall(MatDestroy(&_hc_radial_heat_conduction_mat));
59 102 : LibmeshPetscCall(VecDestroy(&_hc_radial_heat_conduction_rhs));
60 102 : LibmeshPetscCall(MatDestroy(&_hc_sweep_enthalpy_mat));
61 102 : LibmeshPetscCall(VecDestroy(&_hc_sweep_enthalpy_rhs));
62 102 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
63 : }
64 :
65 : void
66 102 : TriSubChannel1PhaseProblem::initializeSolution()
67 : {
68 102 : 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 6 : auto flat_to_flat = _tri_sch_mesh.getFlatToFlat();
73 6 : auto n_rings = _tri_sch_mesh.getNumOfRings();
74 6 : auto pitch = _subchannel_mesh.getPitch();
75 6 : auto pin_diameter = _subchannel_mesh.getPinDiameter();
76 6 : auto wire_diameter = _tri_sch_mesh.getWireDiameter();
77 6 : auto wire_lead_length = _tri_sch_mesh.getWireLeadLength();
78 6 : auto gap = _tri_sch_mesh.getDuctToPinGap();
79 6 : auto z_blockage = _subchannel_mesh.getZBlockage();
80 6 : auto index_blockage = _subchannel_mesh.getIndexBlockage();
81 6 : auto reduction_blockage = _subchannel_mesh.getReductionBlockage();
82 6 : auto theta = std::acos(wire_lead_length /
83 6 : std::sqrt(std::pow(wire_lead_length, 2) +
84 6 : std::pow(libMesh::pi * (pin_diameter + wire_diameter), 2)));
85 312 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
86 : {
87 38862 : 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 312 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
198 : {
199 57222 : 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 6 : }
238 :
239 3210 : 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 102 : _aux->solution().close();
252 102 : }
253 :
254 : Real
255 6405444 : TriSubChannel1PhaseProblem::computeFrictionFactor(FrictionStruct friction_args)
256 : {
257 : // The upgraded Cheng and Todreas correlation for pressure drop in hexagonal wire-wrapped rod
258 : // bundles
259 6405444 : auto Re = friction_args.Re;
260 6405444 : auto i_ch = friction_args.i_ch;
261 6405444 : auto S = friction_args.S;
262 6405444 : auto w_perim = friction_args.w_perim;
263 6405444 : auto Dh_i = 4.0 * S / w_perim;
264 : Real aL, b1L, b2L, cL;
265 : Real aT, b1T, b2T, cT;
266 6405444 : const Real & pitch = _subchannel_mesh.getPitch();
267 6405444 : const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
268 6405444 : const Real & wire_lead_length = _tri_sch_mesh.getWireLeadLength();
269 6405444 : const Real & wire_diameter = _tri_sch_mesh.getWireDiameter();
270 6405444 : auto p_over_d = pitch / pin_diameter;
271 6405444 : 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 6405444 : auto gap = _tri_sch_mesh.getDuctToPinGap();
275 6405444 : auto w_over_d = (pin_diameter + gap) / pin_diameter;
276 6405444 : auto ReL = std::pow(10, (p_over_d - 1)) * 320.0;
277 6405444 : auto ReT = std::pow(10, 0.7 * (p_over_d - 1)) * 1.0E+4;
278 6405444 : auto psi = std::log(Re / ReL) / std::log(ReT / ReL);
279 : const Real lambda = 7.0;
280 6405444 : auto theta = std::acos(wire_lead_length /
281 6405444 : std::sqrt(std::pow(wire_lead_length, 2) +
282 6405444 : std::pow(libMesh::pi * (pin_diameter + wire_diameter), 2)));
283 6405444 : auto wd_t = (19.56 - 98.71 * (wire_diameter / pin_diameter) +
284 6405444 : 303.47 * std::pow((wire_diameter / pin_diameter), 2.0)) *
285 6405444 : std::pow((wire_lead_length / pin_diameter), -0.541);
286 6405444 : auto wd_l = 1.4 * wd_t;
287 6405444 : 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 6405444 : if (subch_type == EChannelType::CENTER)
297 : {
298 3823488 : 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 3823488 : 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 3823488 : cT = aT + b1T * (p_over_d - 1) + b2T * std::pow((p_over_d - 1), 2.0);
320 : }
321 2581956 : else if (subch_type == EChannelType::EDGE)
322 : {
323 1749384 : 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 1749384 : 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 1749384 : cT = aT + b1T * (w_over_d - 1) + b2T * std::pow((w_over_d - 1), 2.0);
345 : }
346 : else
347 : {
348 832572 : 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 832572 : 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 832572 : 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 6405444 : if ((wire_diameter != 0.0) && (wire_lead_length != 0.0))
376 : {
377 6125004 : if (subch_type == EChannelType::CENTER)
378 : {
379 : // wetted perimeter for center subchannel and bare Pin bundle
380 3651408 : pw_p = libMesh::pi * pin_diameter / 2.0;
381 : // wire projected area - center subchannel wire-wrapped bundle
382 3651408 : ar = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
383 : // bare Pin bundle center subchannel flow area (normal area + wire area)
384 3651408 : 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 3651408 : cT *= (pw_p / w_perim);
387 3651408 : cT += wd_t * (3.0 * ar / a_p) * (Dh_i / wire_lead_length) *
388 3651408 : std::pow((Dh_i / wire_diameter), 0.18);
389 : // laminar friction factor equation constant - Center subchannel
390 3651408 : cL *= (pw_p / w_perim);
391 3651408 : cL += wd_l * (3.0 * ar / a_p) * (Dh_i / wire_lead_length) * (Dh_i / wire_diameter);
392 : }
393 2473596 : else if (subch_type == EChannelType::EDGE)
394 : {
395 : // wire projected area - edge subchannel wire-wrapped bundle
396 1675944 : ar = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 4.0;
397 : // bare Pin bundle edge subchannel flow area (normal area + wire area)
398 1675944 : 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 1675944 : 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 1675944 : 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 797652 : ar = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
408 : // bare Pin bundle corner subchannel flow area (normal area + wire area)
409 797652 : 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 797652 : 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 797652 : cL *= (1 + ws_l * (ar / a_p) * std::pow(std::tan(theta), 2.0));
414 : }
415 : }
416 :
417 : // laminar friction factor
418 6405444 : auto fL = cL / Re;
419 : // turbulent friction factor
420 6405444 : auto fT = cT / std::pow(Re, 0.18);
421 :
422 6405444 : if (Re < ReL)
423 : {
424 : // laminar flow
425 : return fL;
426 : }
427 6405444 : 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 16217640 : TriSubChannel1PhaseProblem::computeBeta(unsigned int i_gap, unsigned int iz, bool enthalpy)
442 : {
443 : auto beta = std::numeric_limits<double>::quiet_NaN();
444 16217640 : const Real & pitch = _subchannel_mesh.getPitch();
445 16217640 : const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
446 16217640 : const Real & wire_lead_length = _tri_sch_mesh.getWireLeadLength();
447 16217640 : const Real & wire_diameter = _tri_sch_mesh.getWireDiameter();
448 16217640 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
449 16217640 : auto Nr = _tri_sch_mesh._n_rings;
450 : unsigned int i_ch = chans.first;
451 : unsigned int j_ch = chans.second;
452 16217640 : auto subch_type_i = _subchannel_mesh.getSubchannelType(i_ch);
453 16217640 : auto subch_type_j = _subchannel_mesh.getSubchannelType(j_ch);
454 16217640 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
455 16217640 : auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
456 16217640 : auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
457 16217640 : auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
458 16217640 : auto Si_in = (*_S_flow_soln)(node_in_i);
459 16217640 : auto Sj_in = (*_S_flow_soln)(node_in_j);
460 16217640 : auto Si_out = (*_S_flow_soln)(node_out_i);
461 16217640 : auto Sj_out = (*_S_flow_soln)(node_out_j);
462 16217640 : auto S_total = Si_in + Sj_in + Si_out + Sj_out;
463 16217640 : auto Si = 0.5 * (Si_in + Si_out);
464 16217640 : auto Sj = 0.5 * (Sj_in + Sj_out);
465 16217640 : auto w_perim_i = 0.5 * ((*_w_perim_soln)(node_in_i) + (*_w_perim_soln)(node_out_i));
466 16217640 : auto w_perim_j = 0.5 * ((*_w_perim_soln)(node_in_j) + (*_w_perim_soln)(node_out_j));
467 16217640 : auto avg_mu = (1 / S_total) * ((*_mu_soln)(node_out_i)*Si_out + (*_mu_soln)(node_in_i)*Si_in +
468 16217640 : (*_mu_soln)(node_out_j)*Sj_out + (*_mu_soln)(node_in_j)*Sj_in);
469 16217640 : auto avg_hD = 4.0 * (Si + Sj) / (w_perim_i + w_perim_j);
470 : auto avg_massflux =
471 16217640 : 0.5 * (((*_mdot_soln)(node_in_i) + (*_mdot_soln)(node_in_j)) / (Si_in + Sj_in) +
472 16217640 : ((*_mdot_soln)(node_out_i) + (*_mdot_soln)(node_out_j)) / (Si_out + Sj_out));
473 16217640 : auto Re = avg_massflux * avg_hD / avg_mu;
474 : // crossflow area between channels i,j (dz*gap_width)
475 16217640 : auto gap = _subchannel_mesh.getGapWidth(iz, i_gap);
476 : // Calculation of flow regime
477 16217640 : auto ReL = 320.0 * std::pow(10.0, pitch / pin_diameter - 1);
478 16217640 : 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 16217640 : if ((subch_type_i == EChannelType::CENTER || subch_type_j == EChannelType::CENTER) &&
483 9751320 : (wire_lead_length != 0) && (wire_diameter != 0))
484 : {
485 : // Calculation of geometric parameters
486 : // wire angle
487 9161640 : auto theta = std::acos(wire_lead_length /
488 9161640 : std::sqrt(std::pow(wire_lead_length, 2) +
489 9161640 : std::pow(libMesh::pi * (pin_diameter + wire_diameter), 2)));
490 : // projected area of wire on subchannel
491 9161640 : auto Ar1 = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
492 : // bare subchannel flow area
493 : auto A1prime =
494 9161640 : (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 9161640 : 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 9161640 : 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 9161640 : auto CmT = CmT_constant * std::pow((pitch - pin_diameter) / pin_diameter, -0.5);
514 9161640 : auto CmL = CmL_constant * std::pow((pitch - pin_diameter) / pin_diameter, -0.5);
515 :
516 9161640 : if (Re < ReL)
517 : {
518 : Cm = CmL;
519 : }
520 9161640 : 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 9161640 : beta = Cm * std::pow(Ar1 / A1, 0.5) * std::tan(theta);
532 9161640 : }
533 : // EDGE OR CORNER SUBCHANNELS/ SWEEP FLOW
534 7056000 : else if ((subch_type_i == EChannelType::CORNER || subch_type_i == EChannelType::EDGE) &&
535 6466320 : (subch_type_j == EChannelType::CORNER || subch_type_j == EChannelType::EDGE) &&
536 6466320 : (wire_lead_length != 0) && (wire_diameter != 0))
537 : {
538 6249600 : auto theta = std::acos(wire_lead_length /
539 6249600 : std::sqrt(std::pow(wire_lead_length, 2) +
540 6249600 : std::pow(libMesh::pi * (pin_diameter + wire_diameter), 2)));
541 : // Calculation of geometric parameters
542 : // distance from pin surface to duct
543 6249600 : auto dpgap = _tri_sch_mesh.getDuctToPinGap();
544 : // Edge pitch parameter defined as pin diameter plus distance to duct wall
545 6249600 : auto w = pin_diameter + dpgap;
546 6249600 : auto Ar2 = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 4.0;
547 6249600 : auto A2prime = pitch * (w - pin_diameter / 2.0) - libMesh::pi * std::pow(pin_diameter, 2) / 8.0;
548 6249600 : 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 6249600 : 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 6249600 : auto CsL = CsL_constant * std::pow(wire_lead_length / pin_diameter, 0.3);
564 6249600 : auto CsT = CsT_constant * std::pow(wire_lead_length / pin_diameter, 0.3);
565 :
566 6249600 : if (Re < ReL)
567 : {
568 : Cs = CsL;
569 : }
570 6249600 : 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 6249600 : 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 16217640 : 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 : void
677 2262 : TriSubChannel1PhaseProblem::computeh(int iblock)
678 : {
679 2262 : unsigned int last_node = (iblock + 1) * _block_size;
680 2262 : unsigned int first_node = iblock * _block_size + 1;
681 2262 : const Real & wire_lead_length = _tri_sch_mesh.getWireLeadLength();
682 2262 : const Real & wire_diameter = _tri_sch_mesh.getWireDiameter();
683 2262 : const Real & pitch = _subchannel_mesh.getPitch();
684 2262 : const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
685 :
686 2262 : if (iblock == 0)
687 : {
688 112242 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
689 : {
690 109980 : auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
691 109980 : auto h_out = _fp->h_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
692 109980 : if (h_out < 0)
693 : {
694 0 : mooseError(
695 0 : name(), " : Calculation of negative Enthalpy h_out = : ", h_out, " Axial Level= : ", 0);
696 : }
697 109980 : _h_soln->set(node, h_out);
698 : }
699 : }
700 :
701 2262 : if (!_implicit_bool)
702 : {
703 108 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
704 : {
705 90 : auto z_grid = _subchannel_mesh.getZGrid();
706 90 : auto dz = z_grid[iz] - z_grid[iz - 1];
707 : // Calculation of average mass flux of all periphery subchannels
708 : Real edge_flux_ave = 0.0;
709 : Real mdot_sum = 0.0;
710 : Real si_sum = 0.0;
711 3870 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
712 : {
713 3780 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
714 3780 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
715 : {
716 1620 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
717 1620 : auto Si = (*_S_flow_soln)(node_in);
718 1620 : auto mdot_in = (*_mdot_soln)(node_in);
719 1620 : mdot_sum = mdot_sum + mdot_in;
720 1620 : si_sum = si_sum + Si;
721 : }
722 : }
723 90 : edge_flux_ave = mdot_sum / si_sum;
724 :
725 3870 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
726 : {
727 3780 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
728 3780 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
729 3780 : auto mdot_in = (*_mdot_soln)(node_in);
730 3780 : auto h_in = (*_h_soln)(node_in); // J/kg
731 3780 : auto volume = dz * (*_S_flow_soln)(node_in);
732 3780 : auto mdot_out = (*_mdot_soln)(node_out);
733 3780 : auto h_out = 0.0;
734 : Real sumWijh = 0.0;
735 : Real sumWijPrimeDhij = 0.0;
736 : Real sweep_enthalpy = 0.0;
737 : Real e_cond = 0.0;
738 :
739 : // Calculate added enthalpy from heatflux (Pin, Duct)
740 3780 : Real added_enthalpy = computeAddedHeatPin(i_ch, iz);
741 3780 : added_enthalpy += computeAddedHeatDuct(i_ch, iz);
742 :
743 : // Calculate net sum of enthalpy into/out-of channel i from channels j around i
744 : // (Turbulent diffusion, Diversion Crossflow, Sweep flow Enthalpy, Radial heat conduction)
745 : unsigned int counter = 0;
746 14580 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
747 : {
748 10800 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
749 10800 : auto gap = _subchannel_mesh.getGapWidth(iz, i_gap);
750 10800 : auto Sij = dz * gap;
751 : unsigned int ii_ch = chans.first; // the first subchannel next to gap i_gap
752 : unsigned int jj_ch = chans.second; // the second subchannel next to gap i_gap
753 10800 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
754 10800 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
755 10800 : auto subch_type_i = _subchannel_mesh.getSubchannelType(ii_ch);
756 10800 : auto subch_type_j = _subchannel_mesh.getSubchannelType(jj_ch);
757 : // Define donor enthalpy
758 : auto h_star = 0.0;
759 10800 : if (_Wij(i_gap, iz) > 0.0)
760 8640 : h_star = (*_h_soln)(node_in_i);
761 2160 : else if (_Wij(i_gap, iz) < 0.0)
762 2160 : h_star = (*_h_soln)(node_in_j);
763 : // Diversion crossflow
764 : // take care of the sign by applying the map, use donor cell
765 10800 : sumWijh += _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * h_star;
766 10800 : counter++;
767 : // SWEEP FLOW is calculated if i_gap is located in the periphery
768 : // and we have a wire-wrap (if i_gap is in the periphery then i_ch is in the periphery)
769 : // There are two gaps per periphery subchannel that this is true.
770 10800 : if ((subch_type_i == EChannelType::CORNER || subch_type_i == EChannelType::EDGE) &&
771 3240 : (subch_type_j == EChannelType::CORNER || subch_type_j == EChannelType::EDGE) &&
772 3240 : (wire_lead_length != 0) && (wire_diameter != 0))
773 : {
774 : // donor subchannel and node of sweep flow. The donor subchannel is the subchannel next
775 : // to i_ch that sweep flow, flows from and into i_ch
776 3240 : auto sweep_donor = _tri_sch_mesh.getSweepFlowChans(i_ch).first;
777 3240 : auto * node_sweep_donor = _subchannel_mesh.getChannelNode(sweep_donor, iz - 1);
778 : // if one of the neighbor subchannels of the periphery gap is the donor subchannel
779 : //(the other would be the i_ch) sweep enthalpy flows into i_ch
780 3240 : if ((ii_ch == sweep_donor) || (jj_ch == sweep_donor))
781 : {
782 1620 : sweep_enthalpy +=
783 1620 : computeBeta(i_gap, iz, true) * edge_flux_ave * Sij * (*_h_soln)(node_sweep_donor);
784 : }
785 : // else sweep enthalpy flows out of i_ch
786 : else
787 : {
788 1620 : sweep_enthalpy -=
789 1620 : computeBeta(i_gap, iz, true) * edge_flux_ave * Sij * (*_h_soln)(node_in);
790 : }
791 : }
792 : // Inner gap
793 : // Turbulent Diffusion
794 : else
795 : {
796 7560 : sumWijPrimeDhij +=
797 7560 : _WijPrime(i_gap, iz) * (2 * h_in - (*_h_soln)(node_in_j) - (*_h_soln)(node_in_i));
798 : }
799 :
800 : // compute the radial heat conduction through the gaps
801 : Real dist_ij = pitch;
802 :
803 10800 : if (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::EDGE)
804 : {
805 1080 : dist_ij = pitch;
806 : }
807 9720 : else if ((subch_type_i == EChannelType::CORNER && subch_type_j == EChannelType::EDGE) ||
808 9540 : (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::CORNER))
809 : {
810 2160 : dist_ij = pitch;
811 : }
812 : else
813 : {
814 7560 : dist_ij = pitch / std::sqrt(3);
815 : }
816 :
817 10800 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
818 10800 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
819 : auto shape_factor =
820 10800 : 0.66 * (pitch / pin_diameter) *
821 10800 : std::pow((_subchannel_mesh.getGapWidth(iz, i_gap) / pin_diameter), -0.3);
822 10800 : if (ii_ch == i_ch)
823 : {
824 5400 : e_cond += 0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
825 5400 : ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
826 : }
827 : else
828 : {
829 5400 : e_cond += -0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
830 5400 : ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
831 : }
832 : }
833 :
834 : // compute the axial heat conduction between current and lower axial node
835 3780 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz);
836 3780 : auto * node_in_j = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
837 3780 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
838 3780 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
839 3780 : auto Si = (*_S_flow_soln)(node_in_i);
840 3780 : auto dist_ij = z_grid[iz] - z_grid[iz - 1];
841 :
842 3780 : e_cond += 0.5 * (thcon_i + thcon_j) * Si * ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) /
843 : dist_ij;
844 :
845 3780 : unsigned int nz = _subchannel_mesh.getNumOfAxialCells();
846 : // compute the axial heat conduction between current and upper axial node
847 3780 : if (iz < nz)
848 : {
849 3024 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz);
850 3024 : auto * node_in_j = _subchannel_mesh.getChannelNode(i_ch, iz + 1);
851 3024 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
852 3024 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
853 3024 : auto Si = (*_S_flow_soln)(node_in_i);
854 3024 : auto dist_ij = z_grid[iz + 1] - z_grid[iz];
855 3024 : e_cond += 0.5 * (thcon_i + thcon_j) * Si *
856 3024 : ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
857 : }
858 :
859 : // end of radial heat conduction calc.
860 3780 : h_out =
861 7560 : (mdot_in * h_in - sumWijh - sumWijPrimeDhij + added_enthalpy + e_cond + sweep_enthalpy +
862 3780 : _TR * _rho_soln->old(node_out) * _h_soln->old(node_out) * volume / _dt) /
863 3780 : (mdot_out + _TR * (*_rho_soln)(node_out)*volume / _dt);
864 3780 : if (h_out < 0)
865 : {
866 0 : mooseError(name(),
867 : " : Calculation of negative Enthalpy h_out = : ",
868 : h_out,
869 : " Axial Level= : ",
870 : iz);
871 : }
872 3780 : _h_soln->set(node_out, h_out); // J/kg
873 : }
874 90 : }
875 : }
876 : else
877 : {
878 2244 : LibmeshPetscCall(MatZeroEntries(_hc_time_derivative_mat));
879 2244 : LibmeshPetscCall(MatZeroEntries(_hc_advective_derivative_mat));
880 2244 : LibmeshPetscCall(MatZeroEntries(_hc_cross_derivative_mat));
881 2244 : LibmeshPetscCall(MatZeroEntries(_hc_axial_heat_conduction_mat));
882 2244 : LibmeshPetscCall(MatZeroEntries(_hc_radial_heat_conduction_mat));
883 2244 : LibmeshPetscCall(MatZeroEntries(_hc_sweep_enthalpy_mat));
884 :
885 2244 : LibmeshPetscCall(VecZeroEntries(_hc_time_derivative_rhs));
886 2244 : LibmeshPetscCall(VecZeroEntries(_hc_advective_derivative_rhs));
887 2244 : LibmeshPetscCall(VecZeroEntries(_hc_cross_derivative_rhs));
888 2244 : LibmeshPetscCall(VecZeroEntries(_hc_added_heat_rhs));
889 2244 : LibmeshPetscCall(VecZeroEntries(_hc_axial_heat_conduction_rhs));
890 2244 : LibmeshPetscCall(VecZeroEntries(_hc_radial_heat_conduction_rhs));
891 2244 : LibmeshPetscCall(VecZeroEntries(_hc_sweep_enthalpy_rhs));
892 :
893 2244 : LibmeshPetscCall(MatZeroEntries(_hc_sys_h_mat));
894 2244 : LibmeshPetscCall(VecZeroEntries(_hc_sys_h_rhs));
895 :
896 80844 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
897 : {
898 78600 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
899 78600 : auto pitch = _subchannel_mesh.getPitch();
900 78600 : auto pin_diameter = _subchannel_mesh.getPinDiameter();
901 78600 : auto iz_ind = iz - first_node;
902 :
903 3921240 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
904 : {
905 3842640 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
906 3842640 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
907 3842640 : auto S_in = (*_S_flow_soln)(node_in);
908 3842640 : auto S_out = (*_S_flow_soln)(node_out);
909 3842640 : auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
910 3842640 : auto volume = dz * S_interp;
911 :
912 : PetscScalar Pe = 0.5;
913 3842640 : if (_interpolation_scheme == 3)
914 : {
915 : // Compute the Peclet number
916 942984 : auto w_perim_in = (*_w_perim_soln)(node_in);
917 942984 : auto w_perim_out = (*_w_perim_soln)(node_out);
918 942984 : auto w_perim_interp = this->computeInterpolatedValue(w_perim_out, w_perim_in, 0.5);
919 942984 : auto K_in = _fp->k_from_p_T((*_P_soln)(node_in) + _P_out, (*_T_soln)(node_in));
920 942984 : auto K_out = _fp->k_from_p_T((*_P_soln)(node_out) + _P_out, (*_T_soln)(node_out));
921 942984 : auto K = this->computeInterpolatedValue(K_out, K_in, 0.5);
922 942984 : auto cp_in = _fp->cp_from_p_T((*_P_soln)(node_in) + _P_out, (*_T_soln)(node_in));
923 942984 : auto cp_out = _fp->cp_from_p_T((*_P_soln)(node_out) + _P_out, (*_T_soln)(node_out));
924 942984 : auto cp = this->computeInterpolatedValue(cp_out, cp_in, 0.5);
925 : auto mdot_loc =
926 942984 : this->computeInterpolatedValue((*_mdot_soln)(node_out), (*_mdot_soln)(node_in), 0.5);
927 : // hydraulic diameter in the i direction
928 942984 : auto Dh_i = 4.0 * S_interp / w_perim_interp;
929 942984 : Pe = mdot_loc * Dh_i * cp / (K * S_interp) * (mdot_loc / std::abs(mdot_loc));
930 : }
931 3842640 : auto alpha = computeInterpolationCoefficients(Pe);
932 :
933 : // Time derivative term
934 3842640 : if (iz == first_node)
935 : {
936 : PetscScalar value_vec_tt =
937 109224 : -1.0 * _TR * alpha * (*_rho_soln)(node_in) * (*_h_soln)(node_in)*volume / _dt;
938 109224 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
939 109224 : LibmeshPetscCall(
940 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
941 : }
942 : else
943 : {
944 3733416 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
945 3733416 : PetscInt col_tt = i_ch + _n_channels * (iz_ind - 1);
946 3733416 : PetscScalar value_tt = _TR * alpha * (*_rho_soln)(node_in)*volume / _dt;
947 3733416 : LibmeshPetscCall(MatSetValues(
948 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
949 : }
950 :
951 : // Adding diagonal elements
952 3842640 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
953 3842640 : PetscInt col_tt = i_ch + _n_channels * iz_ind;
954 3842640 : PetscScalar value_tt = _TR * (1.0 - alpha) * (*_rho_soln)(node_out)*volume / _dt;
955 3842640 : LibmeshPetscCall(MatSetValues(
956 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
957 :
958 : // Adding RHS elements
959 : PetscScalar rho_old_interp =
960 3842640 : computeInterpolatedValue(_rho_soln->old(node_out), _rho_soln->old(node_in), Pe);
961 : PetscScalar h_old_interp =
962 3842640 : computeInterpolatedValue(_h_soln->old(node_out), _h_soln->old(node_in), Pe);
963 3842640 : PetscScalar value_vec_tt = _TR * rho_old_interp * h_old_interp * volume / _dt;
964 3842640 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
965 3842640 : LibmeshPetscCall(
966 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
967 :
968 : // Advective derivative term
969 3842640 : if (iz == first_node)
970 : {
971 109224 : PetscInt row_at = i_ch + _n_channels * iz_ind;
972 109224 : PetscScalar value_at = alpha * (*_mdot_soln)(node_in) * (*_h_soln)(node_in);
973 109224 : LibmeshPetscCall(
974 : VecSetValues(_hc_advective_derivative_rhs, 1, &row_at, &value_at, ADD_VALUES));
975 :
976 109224 : value_at = alpha * (*_mdot_soln)(node_out) - (1 - alpha) * (*_mdot_soln)(node_in);
977 109224 : PetscInt col_at = i_ch + _n_channels * iz_ind;
978 109224 : LibmeshPetscCall(MatSetValues(
979 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
980 :
981 109224 : value_at = (1 - alpha) * (*_mdot_soln)(node_out);
982 109224 : col_at = i_ch + _n_channels * (iz_ind + 1);
983 109224 : LibmeshPetscCall(MatSetValues(
984 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
985 : }
986 3733416 : else if (iz == last_node)
987 : {
988 109224 : PetscInt row_at = i_ch + _n_channels * iz_ind;
989 109224 : PetscScalar value_at = 1.0 * (*_mdot_soln)(node_out);
990 109224 : PetscInt col_at = i_ch + _n_channels * iz_ind;
991 109224 : LibmeshPetscCall(MatSetValues(
992 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
993 :
994 109224 : value_at = -1.0 * (*_mdot_soln)(node_in);
995 109224 : col_at = i_ch + _n_channels * (iz_ind - 1);
996 109224 : LibmeshPetscCall(MatSetValues(
997 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
998 : }
999 : else
1000 : {
1001 3624192 : PetscInt row_at = i_ch + _n_channels * iz_ind;
1002 : PetscInt col_at;
1003 :
1004 3624192 : PetscScalar value_at = -alpha * (*_mdot_soln)(node_in);
1005 3624192 : col_at = i_ch + _n_channels * (iz_ind - 1);
1006 3624192 : LibmeshPetscCall(MatSetValues(
1007 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
1008 :
1009 3624192 : value_at = alpha * (*_mdot_soln)(node_out) - (1 - alpha) * (*_mdot_soln)(node_in);
1010 3624192 : col_at = i_ch + _n_channels * iz_ind;
1011 3624192 : LibmeshPetscCall(MatSetValues(
1012 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
1013 :
1014 3624192 : value_at = (1 - alpha) * (*_mdot_soln)(node_out);
1015 3624192 : col_at = i_ch + _n_channels * (iz_ind + 1);
1016 3624192 : LibmeshPetscCall(MatSetValues(
1017 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
1018 : }
1019 :
1020 : // Axial heat conduction
1021 3842640 : auto * node_center = _subchannel_mesh.getChannelNode(i_ch, iz);
1022 3842640 : auto K_center = _fp->k_from_p_T((*_P_soln)(node_center) + _P_out, (*_T_soln)(node_center));
1023 : auto cp_center =
1024 3842640 : _fp->cp_from_p_T((*_P_soln)(node_center) + _P_out, (*_T_soln)(node_center));
1025 3842640 : auto diff_center = K_center / (cp_center + 1e-15);
1026 :
1027 3842640 : if (iz == first_node)
1028 : {
1029 109224 : auto * node_top = _subchannel_mesh.getChannelNode(i_ch, iz + 1);
1030 109224 : auto * node_bottom = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1031 : auto K_bottom =
1032 109224 : _fp->k_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1033 109224 : auto K_top = _fp->k_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
1034 : auto cp_bottom =
1035 109224 : _fp->cp_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1036 109224 : auto cp_top = _fp->cp_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
1037 109224 : auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
1038 109224 : auto diff_top = K_top / (cp_top + 1e-15);
1039 :
1040 109224 : auto dz_up = _z_grid[iz + 1] - _z_grid[iz];
1041 109224 : auto dz_down = _z_grid[iz] - _z_grid[iz - 1];
1042 : auto S_up =
1043 109224 : computeInterpolatedValue((*_S_flow_soln)(node_top), (*_S_flow_soln)(node_center));
1044 : auto S_down =
1045 109224 : computeInterpolatedValue((*_S_flow_soln)(node_center), (*_S_flow_soln)(node_bottom));
1046 109224 : auto diff_up = computeInterpolatedValue(diff_top, diff_center);
1047 109224 : auto diff_down = computeInterpolatedValue(diff_center, diff_bottom);
1048 :
1049 : // Diagonal value
1050 109224 : PetscInt row_at = i_ch + _n_channels * iz_ind;
1051 109224 : PetscInt col_at = i_ch + _n_channels * iz_ind;
1052 109224 : PetscScalar value_at = diff_up * S_up / dz_up + diff_down * S_down / dz_down;
1053 109224 : LibmeshPetscCall(MatSetValues(
1054 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1055 :
1056 : // Bottom value
1057 109224 : value_at = 1.0 * diff_down * S_down / dz_down * (*_h_soln)(node_bottom);
1058 109224 : LibmeshPetscCall(
1059 : VecSetValues(_hc_axial_heat_conduction_rhs, 1, &row_at, &value_at, ADD_VALUES));
1060 :
1061 : // Top value
1062 109224 : col_at = i_ch + _n_channels * (iz_ind + 1);
1063 109224 : value_at = -diff_up * S_up / dz_up;
1064 109224 : LibmeshPetscCall(MatSetValues(
1065 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1066 : }
1067 3733416 : else if (iz == last_node)
1068 : {
1069 109224 : auto * node_bottom = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1070 : auto K_bottom =
1071 109224 : _fp->k_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1072 : auto cp_bottom =
1073 109224 : _fp->cp_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1074 109224 : auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
1075 :
1076 109224 : auto dz_down = _z_grid[iz] - _z_grid[iz - 1];
1077 109224 : auto S_down = 0.5 * ((*_S_flow_soln)(node_center) + (*_S_flow_soln)(node_bottom));
1078 109224 : auto diff_down = 0.5 * (diff_center + diff_bottom);
1079 :
1080 : // Diagonal value
1081 109224 : PetscInt row_at = i_ch + _n_channels * iz_ind;
1082 109224 : PetscInt col_at = i_ch + _n_channels * iz_ind;
1083 109224 : PetscScalar value_at = diff_down * S_down / dz_down;
1084 109224 : LibmeshPetscCall(MatSetValues(
1085 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1086 :
1087 : // Bottom value
1088 109224 : col_at = i_ch + _n_channels * (iz_ind - 1);
1089 109224 : value_at = -diff_down * S_down / dz_down;
1090 109224 : LibmeshPetscCall(MatSetValues(
1091 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1092 :
1093 : // Outflow derivative
1094 : /// TODO: Current axial derivative is zero - check if outflow conditions may make a difference
1095 : // value_at = -1.0 * (*_mdot_soln)(node_center) * (*_h_soln)(node_center);
1096 : // VecSetValues(_hc_axial_heat_conduction_rhs, 1, &row_at, &value_at, ADD_VALUES);
1097 : }
1098 : else
1099 : {
1100 3624192 : auto * node_top = _subchannel_mesh.getChannelNode(i_ch, iz + 1);
1101 3624192 : auto * node_bottom = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1102 : auto K_bottom =
1103 3624192 : _fp->k_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1104 3624192 : auto K_top = _fp->k_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
1105 : auto cp_bottom =
1106 3624192 : _fp->cp_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
1107 3624192 : auto cp_top = _fp->cp_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
1108 3624192 : auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
1109 3624192 : auto diff_top = K_top / (cp_top + 1e-15);
1110 :
1111 3624192 : auto dz_up = _z_grid[iz + 1] - _z_grid[iz];
1112 3624192 : auto dz_down = _z_grid[iz] - _z_grid[iz - 1];
1113 : auto S_up =
1114 3624192 : computeInterpolatedValue((*_S_flow_soln)(node_top), (*_S_flow_soln)(node_center));
1115 : auto S_down =
1116 3624192 : computeInterpolatedValue((*_S_flow_soln)(node_center), (*_S_flow_soln)(node_bottom));
1117 3624192 : auto diff_up = computeInterpolatedValue(diff_top, diff_center);
1118 3624192 : auto diff_down = computeInterpolatedValue(diff_center, diff_bottom);
1119 :
1120 : // Diagonal value
1121 3624192 : PetscInt row_at = i_ch + _n_channels * iz_ind;
1122 3624192 : PetscInt col_at = i_ch + _n_channels * iz_ind;
1123 3624192 : PetscScalar value_at = diff_up * S_up / dz_up + diff_down * S_down / dz_down;
1124 3624192 : LibmeshPetscCall(MatSetValues(
1125 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1126 :
1127 : // Bottom value
1128 3624192 : col_at = i_ch + _n_channels * (iz_ind - 1);
1129 3624192 : value_at = -diff_down * S_down / dz_down;
1130 3624192 : LibmeshPetscCall(MatSetValues(
1131 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1132 :
1133 : // Top value
1134 3624192 : col_at = i_ch + _n_channels * (iz_ind + 1);
1135 3624192 : value_at = -diff_up * S_up / dz_up;
1136 3624192 : LibmeshPetscCall(MatSetValues(
1137 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
1138 : }
1139 :
1140 : // Radial Terms
1141 : unsigned int counter = 0;
1142 : unsigned int cross_index = iz;
1143 : // Real radial_heat_conduction(0.0);
1144 14898960 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
1145 : {
1146 11056320 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
1147 : unsigned int ii_ch = chans.first;
1148 : unsigned int jj_ch = chans.second;
1149 11056320 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
1150 11056320 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
1151 : PetscScalar h_star;
1152 : // figure out donor axial velocity
1153 11056320 : if (_Wij(i_gap, cross_index) > 0.0)
1154 : {
1155 6766044 : if (iz == first_node)
1156 : {
1157 208644 : h_star = (*_h_soln)(node_in_i);
1158 417288 : PetscScalar value_vec_ct = -1.0 * alpha *
1159 208644 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1160 208644 : _Wij(i_gap, cross_index) * h_star;
1161 208644 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
1162 208644 : LibmeshPetscCall(VecSetValues(
1163 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
1164 : }
1165 : else
1166 : {
1167 6557400 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1168 6557400 : _Wij(i_gap, cross_index);
1169 6557400 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1170 6557400 : PetscInt col_ct = ii_ch + _n_channels * (iz_ind - 1);
1171 6557400 : LibmeshPetscCall(MatSetValues(
1172 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
1173 : }
1174 6766044 : PetscScalar value_ct = (1.0 - alpha) *
1175 6766044 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1176 6766044 : _Wij(i_gap, cross_index);
1177 6766044 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1178 6766044 : PetscInt col_ct = ii_ch + _n_channels * iz_ind;
1179 6766044 : LibmeshPetscCall(MatSetValues(
1180 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
1181 : }
1182 4290276 : else if (_Wij(i_gap, cross_index) < 0.0) // _Wij=0 operations not necessary
1183 : {
1184 3941184 : if (iz == first_node)
1185 : {
1186 103968 : h_star = (*_h_soln)(node_in_j);
1187 207936 : PetscScalar value_vec_ct = -1.0 * alpha *
1188 103968 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1189 103968 : _Wij(i_gap, cross_index) * h_star;
1190 103968 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
1191 103968 : LibmeshPetscCall(VecSetValues(
1192 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
1193 : }
1194 : else
1195 : {
1196 3837216 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1197 3837216 : _Wij(i_gap, cross_index);
1198 3837216 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1199 3837216 : PetscInt col_ct = jj_ch + _n_channels * (iz_ind - 1);
1200 3837216 : LibmeshPetscCall(MatSetValues(
1201 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
1202 : }
1203 3941184 : PetscScalar value_ct = (1.0 - alpha) *
1204 3941184 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
1205 3941184 : _Wij(i_gap, cross_index);
1206 3941184 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1207 3941184 : PetscInt col_ct = jj_ch + _n_channels * iz_ind;
1208 3941184 : LibmeshPetscCall(MatSetValues(
1209 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
1210 : }
1211 :
1212 : // Turbulent cross flows
1213 11056320 : if (iz == first_node)
1214 : {
1215 : PetscScalar value_vec_ct =
1216 314208 : -2.0 * alpha * (*_h_soln)(node_in)*_WijPrime(i_gap, cross_index);
1217 314208 : value_vec_ct += alpha * (*_h_soln)(node_in_j)*_WijPrime(i_gap, cross_index);
1218 314208 : value_vec_ct += alpha * (*_h_soln)(node_in_i)*_WijPrime(i_gap, cross_index);
1219 314208 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
1220 314208 : LibmeshPetscCall(
1221 : VecSetValues(_hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
1222 : }
1223 : else
1224 : {
1225 10742112 : PetscScalar value_center_ct = 2.0 * alpha * _WijPrime(i_gap, cross_index);
1226 10742112 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1227 10742112 : PetscInt col_ct = i_ch + _n_channels * (iz_ind - 1);
1228 10742112 : LibmeshPetscCall(MatSetValues(
1229 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
1230 :
1231 10742112 : PetscScalar value_left_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
1232 10742112 : row_ct = i_ch + _n_channels * iz_ind;
1233 10742112 : col_ct = jj_ch + _n_channels * (iz_ind - 1);
1234 10742112 : LibmeshPetscCall(MatSetValues(
1235 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
1236 :
1237 10742112 : PetscScalar value_right_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
1238 10742112 : row_ct = i_ch + _n_channels * iz_ind;
1239 10742112 : col_ct = ii_ch + _n_channels * (iz_ind - 1);
1240 10742112 : LibmeshPetscCall(MatSetValues(
1241 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
1242 : }
1243 11056320 : PetscScalar value_center_ct = 2.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
1244 11056320 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
1245 11056320 : PetscInt col_ct = i_ch + _n_channels * iz_ind;
1246 11056320 : LibmeshPetscCall(MatSetValues(
1247 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
1248 :
1249 11056320 : PetscScalar value_left_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
1250 11056320 : row_ct = i_ch + _n_channels * iz_ind;
1251 11056320 : col_ct = jj_ch + _n_channels * iz_ind;
1252 11056320 : LibmeshPetscCall(MatSetValues(
1253 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
1254 :
1255 11056320 : PetscScalar value_right_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
1256 11056320 : row_ct = i_ch + _n_channels * iz_ind;
1257 11056320 : col_ct = ii_ch + _n_channels * iz_ind;
1258 11056320 : LibmeshPetscCall(MatSetValues(
1259 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
1260 :
1261 : // Radial heat conduction
1262 11056320 : auto subch_type_i = _subchannel_mesh.getSubchannelType(ii_ch);
1263 11056320 : auto subch_type_j = _subchannel_mesh.getSubchannelType(jj_ch);
1264 : Real dist_ij = pitch;
1265 :
1266 11056320 : if (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::EDGE)
1267 : {
1268 : dist_ij = pitch;
1269 : }
1270 9951840 : else if ((subch_type_i == EChannelType::CORNER && subch_type_j == EChannelType::EDGE) ||
1271 9794640 : (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::CORNER))
1272 : {
1273 : dist_ij = pitch;
1274 : }
1275 : else
1276 : {
1277 8065440 : dist_ij = pitch / std::sqrt(3);
1278 : }
1279 :
1280 11056320 : auto Sij = dz * _subchannel_mesh.getGapWidth(iz, i_gap);
1281 11056320 : auto K_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
1282 11056320 : auto K_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
1283 11056320 : auto cp_i = _fp->cp_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
1284 11056320 : auto cp_j = _fp->cp_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
1285 11056320 : auto A_i = K_i / cp_i;
1286 11056320 : auto A_j = K_j / cp_j;
1287 11056320 : auto harm_A = 2.0 * A_i * A_j / (A_i + A_j);
1288 : auto shape_factor =
1289 11056320 : 0.66 * (pitch / pin_diameter) *
1290 11056320 : std::pow((_subchannel_mesh.getGapWidth(iz, i_gap) / pin_diameter), -0.3);
1291 : // auto base_value = 0.5 * (A_i + A_j) * Sij * shape_factor / dist_ij;
1292 11056320 : auto base_value = harm_A * shape_factor * Sij / dist_ij;
1293 11056320 : auto neg_base_value = -1.0 * base_value;
1294 :
1295 11056320 : row_ct = ii_ch + _n_channels * iz_ind;
1296 11056320 : col_ct = ii_ch + _n_channels * iz_ind;
1297 11056320 : LibmeshPetscCall(MatSetValues(
1298 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &base_value, ADD_VALUES));
1299 :
1300 11056320 : row_ct = jj_ch + _n_channels * iz_ind;
1301 11056320 : col_ct = jj_ch + _n_channels * iz_ind;
1302 11056320 : LibmeshPetscCall(MatSetValues(
1303 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &base_value, ADD_VALUES));
1304 :
1305 11056320 : row_ct = ii_ch + _n_channels * iz_ind;
1306 11056320 : col_ct = jj_ch + _n_channels * iz_ind;
1307 11056320 : LibmeshPetscCall(MatSetValues(
1308 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &neg_base_value, ADD_VALUES));
1309 :
1310 11056320 : row_ct = jj_ch + _n_channels * iz_ind;
1311 11056320 : col_ct = ii_ch + _n_channels * iz_ind;
1312 11056320 : LibmeshPetscCall(MatSetValues(
1313 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &neg_base_value, ADD_VALUES));
1314 11056320 : counter++;
1315 : }
1316 :
1317 : // Compute the sweep flow enthalpy change
1318 : // Calculation of average mass flux of all periphery subchannels
1319 : Real edge_flux_ave = 0.0;
1320 : Real mdot_sum = 0.0;
1321 : Real si_sum = 0.0;
1322 226819440 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1323 : {
1324 222976800 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
1325 222976800 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
1326 : {
1327 78222240 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1328 78222240 : auto Si = (*_S_flow_soln)(node_in);
1329 78222240 : auto mdot_in = (*_mdot_soln)(node_in);
1330 78222240 : mdot_sum = mdot_sum + mdot_in;
1331 78222240 : si_sum = si_sum + Si;
1332 : }
1333 : }
1334 3842640 : edge_flux_ave = mdot_sum / si_sum;
1335 3842640 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
1336 : PetscScalar sweep_enthalpy = 0.0;
1337 3842640 : if ((subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER) &&
1338 1495440 : (wire_diameter != 0.0) && (wire_lead_length != 0.0))
1339 : {
1340 : auto beta_in = std::numeric_limits<double>::quiet_NaN();
1341 : auto beta_out = std::numeric_limits<double>::quiet_NaN();
1342 : // donor sweep channel for i_ch
1343 1392840 : auto sweep_donor = _tri_sch_mesh.getSweepFlowChans(i_ch).first;
1344 1392840 : auto * node_sweep_donor = _subchannel_mesh.getChannelNode(sweep_donor, iz - 1);
1345 : // Calculation of turbulent mixing parameter
1346 5133960 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
1347 : {
1348 3741120 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
1349 : unsigned int ii_ch = chans.first;
1350 : unsigned int jj_ch = chans.second;
1351 3741120 : auto subch_type_i = _subchannel_mesh.getSubchannelType(ii_ch);
1352 3741120 : auto subch_type_j = _subchannel_mesh.getSubchannelType(jj_ch);
1353 3741120 : if ((subch_type_i == EChannelType::CORNER || subch_type_i == EChannelType::EDGE) &&
1354 2785680 : (subch_type_j == EChannelType::CORNER || subch_type_j == EChannelType::EDGE))
1355 : {
1356 2785680 : if ((ii_ch == sweep_donor) || (jj_ch == sweep_donor))
1357 : {
1358 1392840 : beta_in = computeBeta(i_gap, iz, true);
1359 : }
1360 : else
1361 : {
1362 1392840 : beta_out = computeBeta(i_gap, iz, true);
1363 : }
1364 : }
1365 : }
1366 : // Abort execution if required values are unset
1367 : mooseAssert(!std::isnan(beta_in),
1368 : "beta_in was not set. Check gap logic for i_ch = " + std::to_string(i_ch) +
1369 : ", iz = " + std::to_string(iz));
1370 : mooseAssert(!std::isnan(beta_out),
1371 : "beta_out was not set. Check gap logic for i_ch = " + std::to_string(i_ch) +
1372 : ", iz = " + std::to_string(iz));
1373 :
1374 1392840 : auto gap = _tri_sch_mesh.getDuctToPinGap();
1375 1392840 : auto Sij = dz * gap;
1376 1392840 : auto wsweep_in = edge_flux_ave * beta_in * Sij;
1377 1392840 : auto wsweep_out = edge_flux_ave * beta_out * Sij;
1378 1392840 : auto sweep_hin = (*_h_soln)(node_sweep_donor);
1379 1392840 : auto sweep_hout = (*_h_soln)(node_in);
1380 1392840 : sweep_enthalpy = (wsweep_in * sweep_hin - wsweep_out * sweep_hout);
1381 :
1382 1392840 : if (iz == first_node)
1383 : {
1384 40644 : PetscInt row_sh = i_ch + _n_channels * iz_ind;
1385 40644 : PetscScalar value_hs = -sweep_enthalpy;
1386 40644 : LibmeshPetscCall(
1387 : VecSetValues(_hc_sweep_enthalpy_rhs, 1, &row_sh, &value_hs, ADD_VALUES));
1388 : }
1389 : else
1390 : {
1391 : // coefficient of sweep_hin
1392 1352196 : PetscInt row_sh = i_ch + _n_channels * (iz_ind - 1);
1393 1352196 : PetscInt col_sh = i_ch + _n_channels * (iz_ind - 1);
1394 1352196 : LibmeshPetscCall(MatSetValues(
1395 : _hc_sweep_enthalpy_mat, 1, &row_sh, 1, &col_sh, &wsweep_out, ADD_VALUES));
1396 1352196 : PetscInt col_sh_l = sweep_donor + _n_channels * (iz_ind - 1);
1397 1352196 : PetscScalar neg_sweep_in = -1.0 * wsweep_in;
1398 : // coefficient of sweep_hout
1399 1352196 : LibmeshPetscCall(MatSetValues(
1400 : _hc_sweep_enthalpy_mat, 1, &row_sh, 1, &col_sh_l, &(neg_sweep_in), ADD_VALUES));
1401 : }
1402 : }
1403 :
1404 : // Add heat enthalpy from pin and/or duct
1405 3842640 : PetscScalar added_enthalpy = computeAddedHeatPin(i_ch, iz);
1406 3842640 : added_enthalpy += computeAddedHeatDuct(i_ch, iz);
1407 3842640 : PetscInt row_vec_ht = i_ch + _n_channels * iz_ind;
1408 3842640 : LibmeshPetscCall(
1409 : VecSetValues(_hc_added_heat_rhs, 1, &row_vec_ht, &added_enthalpy, ADD_VALUES));
1410 : }
1411 : }
1412 : // Assembling system
1413 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
1414 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
1415 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
1416 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
1417 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
1418 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
1419 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_axial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1420 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_axial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1421 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_radial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1422 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_radial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1423 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sweep_enthalpy_mat, MAT_FINAL_ASSEMBLY));
1424 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sweep_enthalpy_mat, MAT_FINAL_ASSEMBLY));
1425 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1426 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1427 : // Add all matrices together
1428 2244 : LibmeshPetscCall(
1429 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1430 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1431 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1432 2244 : LibmeshPetscCall(
1433 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1434 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1435 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1436 2244 : LibmeshPetscCall(
1437 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1438 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1439 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1440 2244 : LibmeshPetscCall(
1441 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_axial_heat_conduction_mat, DIFFERENT_NONZERO_PATTERN));
1442 2244 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1443 2244 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1444 2244 : LibmeshPetscCall(
1445 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_radial_heat_conduction_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_sweep_enthalpy_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 : if (_verbose_subchannel)
1453 1740 : _console << "Block: " << iblock << " - Enthalpy conservation matrix assembled" << std::endl;
1454 : // RHS
1455 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_time_derivative_rhs));
1456 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_advective_derivative_rhs));
1457 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_cross_derivative_rhs));
1458 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_added_heat_rhs));
1459 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_axial_heat_conduction_rhs));
1460 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_radial_heat_conduction_rhs));
1461 2244 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_sweep_enthalpy_rhs));
1462 :
1463 2244 : if (_segregated_bool || (!_monolithic_thermal_bool))
1464 : {
1465 : // Assembly the matrix system
1466 : KSP ksploc;
1467 : PC pc;
1468 : Vec sol;
1469 1380 : LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &sol));
1470 1380 : LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
1471 1380 : LibmeshPetscCall(KSPSetOperators(ksploc, _hc_sys_h_mat, _hc_sys_h_mat));
1472 1380 : LibmeshPetscCall(KSPGetPC(ksploc, &pc));
1473 1380 : LibmeshPetscCall(PCSetType(pc, PCJACOBI));
1474 1380 : LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit));
1475 1380 : LibmeshPetscCall(KSPSetOptionsPrefix(ksploc, "h_sys_"));
1476 1380 : LibmeshPetscCall(KSPSetFromOptions(ksploc));
1477 1380 : LibmeshPetscCall(KSPSolve(ksploc, _hc_sys_h_rhs, sol));
1478 : // VecView(sol, PETSC_VIEWER_STDOUT_WORLD);
1479 : PetscScalar * xx;
1480 1380 : LibmeshPetscCall(VecGetArray(sol, &xx));
1481 51300 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1482 : {
1483 49920 : auto iz_ind = iz - first_node;
1484 2688000 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1485 : {
1486 2638080 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
1487 2638080 : auto h_out = xx[iz_ind * _n_channels + i_ch];
1488 2638080 : if (h_out < 0)
1489 : {
1490 0 : mooseError(name(),
1491 : " : Calculation of negative Enthalpy h_out = : ",
1492 : h_out,
1493 : " Axial Level= : ",
1494 : iz);
1495 : }
1496 2638080 : _h_soln->set(node_out, h_out);
1497 : }
1498 : }
1499 1380 : LibmeshPetscCall(KSPDestroy(&ksploc));
1500 1380 : LibmeshPetscCall(VecDestroy(&sol));
1501 : }
1502 : }
1503 2262 : }
|
