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 "SubChannel1PhaseProblem.h"
14 : #include "SinglePhaseFluidProperties.h"
15 : #include "SCM.h"
16 : #include <limits> // for std::numeric_limits
17 : #include <cmath> // for std::isnan
18 : #include "SCMMixingClosureBase.h"
19 :
20 : registerMooseObject("SubChannelApp", TriSubChannel1PhaseProblem);
21 :
22 : InputParameters
23 310 : TriSubChannel1PhaseProblem::validParams()
24 : {
25 310 : InputParameters params = SubChannel1PhaseProblem::validParams();
26 310 : params.addClassDescription("Solver class for subchannels in a triangular lattice assembly and "
27 : "bare/wire-wrapped fuel pins");
28 310 : return params;
29 0 : }
30 :
31 155 : TriSubChannel1PhaseProblem::TriSubChannel1PhaseProblem(const InputParameters & params)
32 : : SubChannel1PhaseProblem(params),
33 155 : _tri_sch_mesh(SCM::getMesh<TriSubChannelMesh>(_subchannel_mesh))
34 : {
35 : // Initializing heat conduction system
36 155 : LibmeshPetscCall(createPetscMatrix(
37 : _hc_axial_heat_conduction_mat, _block_size * _n_channels, _block_size * _n_channels));
38 155 : LibmeshPetscCall(createPetscVector(_hc_axial_heat_conduction_rhs, _block_size * _n_channels));
39 155 : LibmeshPetscCall(createPetscMatrix(
40 : _hc_radial_heat_conduction_mat, _block_size * _n_channels, _block_size * _n_channels));
41 155 : LibmeshPetscCall(createPetscVector(_hc_radial_heat_conduction_rhs, _block_size * _n_channels));
42 155 : LibmeshPetscCall(createPetscMatrix(
43 : _hc_sweep_enthalpy_mat, _block_size * _n_channels, _block_size * _n_channels));
44 155 : LibmeshPetscCall(createPetscVector(_hc_sweep_enthalpy_rhs, _block_size * _n_channels));
45 155 : }
46 :
47 152 : TriSubChannel1PhaseProblem::~TriSubChannel1PhaseProblem()
48 : {
49 152 : PetscErrorCode ierr = cleanUp();
50 152 : if (ierr)
51 0 : mooseError(name(), ": Error in memory cleanup");
52 152 : }
53 :
54 : PetscErrorCode
55 152 : TriSubChannel1PhaseProblem::cleanUp()
56 : {
57 : PetscFunctionBegin;
58 : // Clean up heat conduction system
59 152 : LibmeshPetscCall(MatDestroy(&_hc_axial_heat_conduction_mat));
60 152 : LibmeshPetscCall(VecDestroy(&_hc_axial_heat_conduction_rhs));
61 152 : LibmeshPetscCall(MatDestroy(&_hc_radial_heat_conduction_mat));
62 152 : LibmeshPetscCall(VecDestroy(&_hc_radial_heat_conduction_rhs));
63 152 : LibmeshPetscCall(MatDestroy(&_hc_sweep_enthalpy_mat));
64 152 : LibmeshPetscCall(VecDestroy(&_hc_sweep_enthalpy_rhs));
65 152 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
66 : }
67 :
68 : void
69 151 : TriSubChannel1PhaseProblem::initializeSolution()
70 : {
71 151 : detectDeformation();
72 :
73 151 : if (_deformation)
74 : {
75 : // update surface area, wetted perimeter based on: Dpin, displacement
76 : Real standard_area, wire_area, additional_area, wetted_perimeter, displaced_area;
77 5 : auto flat_to_flat = _tri_sch_mesh.getFlatToFlat();
78 5 : auto n_rings = _tri_sch_mesh.getNumOfRings();
79 5 : auto pitch = _subchannel_mesh.getPitch();
80 5 : auto pin_diameter = _subchannel_mesh.getPinDiameter();
81 5 : auto wire_diameter = _tri_sch_mesh.getWireDiameter();
82 5 : auto wire_lead_length = _tri_sch_mesh.getWireLeadLength();
83 5 : auto gap = _tri_sch_mesh.getDuctToPinGap();
84 5 : auto z_blockage = _subchannel_mesh.getZBlockage();
85 5 : auto index_blockage = _subchannel_mesh.getIndexBlockage();
86 5 : auto reduction_blockage = _subchannel_mesh.getReductionBlockage();
87 : auto theta =
88 5 : std::acos(wire_lead_length /
89 5 : std::sqrt(Utility::pow<2>(wire_lead_length) +
90 5 : Utility::pow<2>(libMesh::pi * (pin_diameter + wire_diameter))));
91 35 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
92 : {
93 1290 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
94 : {
95 1260 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
96 1260 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
97 1260 : auto Z = _z_grid[iz];
98 : Real rod_area = 0.0;
99 : Real rod_perimeter = 0.0;
100 4320 : for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
101 : {
102 3060 : auto * pin_node = _subchannel_mesh.getPinNode(i_pin, iz);
103 3060 : if (subch_type == EChannelType::CENTER || subch_type == EChannelType::CORNER)
104 : {
105 2340 : rod_area +=
106 2340 : (1.0 / 6.0) * 0.25 * M_PI * (*_Dpin_soln)(pin_node) * (*_Dpin_soln)(pin_node);
107 2340 : rod_perimeter += (1.0 / 6.0) * M_PI * (*_Dpin_soln)(pin_node);
108 : }
109 : else
110 : {
111 720 : rod_area +=
112 720 : (1.0 / 4.0) * 0.25 * M_PI * (*_Dpin_soln)(pin_node) * (*_Dpin_soln)(pin_node);
113 720 : rod_perimeter += (1.0 / 4.0) * M_PI * (*_Dpin_soln)(pin_node);
114 : }
115 : }
116 :
117 1260 : if (subch_type == EChannelType::CENTER)
118 : {
119 720 : standard_area = Utility::pow<2>(pitch) * std::sqrt(3.0) / 4.0;
120 : additional_area = 0.0;
121 : displaced_area = 0.0;
122 720 : wire_area = libMesh::pi * Utility::pow<2>(wire_diameter) / 8.0 / std::cos(theta);
123 720 : wetted_perimeter = rod_perimeter + 0.5 * libMesh::pi * wire_diameter / std::cos(theta);
124 : }
125 540 : else if (subch_type == EChannelType::EDGE)
126 : {
127 360 : standard_area = pitch * (pin_diameter / 2.0 + gap);
128 : additional_area = 0.0;
129 360 : displaced_area = (*_displacement_soln)(node)*pitch;
130 360 : wire_area = libMesh::pi * Utility::pow<2>(wire_diameter) / 8.0 / std::cos(theta);
131 360 : wetted_perimeter =
132 360 : rod_perimeter + 0.5 * libMesh::pi * wire_diameter / std::cos(theta) + pitch;
133 : }
134 : else
135 : {
136 180 : standard_area = 1.0 / std::sqrt(3.0) * Utility::pow<2>(pin_diameter / 2.0 + gap);
137 : additional_area = 0.0;
138 180 : displaced_area = 1.0 / std::sqrt(3.0) *
139 180 : (pin_diameter + 2.0 * gap + (*_displacement_soln)(node)) *
140 180 : (*_displacement_soln)(node);
141 180 : wire_area = libMesh::pi / 24.0 * Utility::pow<2>(wire_diameter) / std::cos(theta);
142 180 : wetted_perimeter =
143 180 : rod_perimeter + libMesh::pi * wire_diameter / std::cos(theta) / 6.0 +
144 180 : 2.0 / std::sqrt(3.0) * (pin_diameter / 2.0 + gap + (*_displacement_soln)(node));
145 : }
146 :
147 : // Calculate subchannel area
148 1260 : auto subchannel_area =
149 1260 : standard_area + additional_area + displaced_area - rod_area - wire_area;
150 :
151 : // Correct subchannel area and wetted perimeter in case of overlapping pins
152 : auto overlapping_pin_area = 0.0;
153 : auto overlapping_wetted_perimeter = 0.0;
154 4860 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
155 : {
156 3600 : auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
157 : auto pin_1 = gap_pins.first;
158 : auto pin_2 = gap_pins.second;
159 3600 : auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
160 3600 : auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
161 3600 : auto Diameter1 = (*_Dpin_soln)(pin_node_1);
162 3600 : auto Radius1 = Diameter1 / 2.0;
163 3600 : auto Diameter2 = (*_Dpin_soln)(pin_node_2);
164 3600 : auto Radius2 = Diameter2 / 2.0;
165 3600 : auto pitch = _subchannel_mesh.getPitch();
166 :
167 3600 : if (pitch < (Radius1 + Radius2)) // overlapping pins
168 : {
169 0 : mooseWarning(" The gap of index : '", i_gap, " at axial cell ", iz, " ' is blocked.");
170 0 : auto cos1 =
171 0 : (pitch * pitch + Radius1 * Radius1 - Radius2 * Radius2) / (2 * pitch * Radius1);
172 0 : auto cos2 =
173 0 : (pitch * pitch + Radius2 * Radius2 - Radius1 * Radius1) / (2 * pitch * Radius2);
174 0 : auto angle1 = 2.0 * acos(cos1);
175 0 : auto angle2 = 2.0 * acos(cos2);
176 : // half of the intersecting arc-length
177 0 : overlapping_wetted_perimeter += 0.5 * angle1 * Radius1 + 0.5 * angle2 * Radius2;
178 : // Half of the overlapping area
179 0 : overlapping_pin_area +=
180 0 : 0.5 * Radius1 * Radius1 * acos(cos1) + 0.5 * Radius2 * Radius2 * acos(cos2) -
181 0 : 0.25 * sqrt((-pitch + Radius1 + Radius2) * (pitch + Radius1 - Radius2) *
182 0 : (pitch - Radius1 + Radius2) * (pitch + Radius1 + Radius2));
183 : }
184 : }
185 1260 : subchannel_area += overlapping_pin_area; // correct surface area
186 1260 : wetted_perimeter += -overlapping_wetted_perimeter; // correct wetted perimeter
187 :
188 : // Apply area reduction on subchannels affected by blockage
189 : auto index = 0;
190 2520 : for (const auto & i_blockage : index_blockage)
191 : {
192 1260 : if (i_ch == i_blockage && (Z >= z_blockage.front() && Z <= z_blockage.back()))
193 : {
194 5 : subchannel_area *= reduction_blockage[index];
195 : }
196 1260 : index++;
197 : }
198 1260 : _S_flow_soln->set(node, subchannel_area);
199 : _w_perim_soln->set(node, wetted_perimeter);
200 : }
201 : }
202 : // update map of gap between pins (gij) based on: Dpin, displacement
203 35 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
204 : {
205 1830 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
206 : {
207 1800 : auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
208 : auto pin_1 = gap_pins.first;
209 : auto pin_2 = gap_pins.second;
210 1800 : auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
211 1800 : auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
212 :
213 1800 : if (pin_1 == pin_2) // Corner or edge gap
214 : {
215 : auto displacement = 0.0;
216 : auto counter = 0.0;
217 3240 : for (auto i_ch : _subchannel_mesh.getPinChannels(pin_1))
218 : {
219 2700 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
220 2700 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
221 2700 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
222 : {
223 1440 : displacement += (*_displacement_soln)(node);
224 1440 : counter += 1.0;
225 : }
226 : }
227 540 : displacement = displacement / counter;
228 540 : _tri_sch_mesh._gij_map[iz][i_gap] =
229 540 : 0.5 * (flat_to_flat - (n_rings - 1) * pitch * std::sqrt(3.0) -
230 540 : (*_Dpin_soln)(pin_node_1)) +
231 : displacement;
232 : }
233 : else // center gap
234 : {
235 1260 : _tri_sch_mesh._gij_map[iz][i_gap] =
236 1260 : pitch - (*_Dpin_soln)(pin_node_1) / 2.0 - (*_Dpin_soln)(pin_node_2) / 2.0;
237 : }
238 : // if pins come in contact, the gap is zero
239 1800 : if (_tri_sch_mesh._gij_map[iz][i_gap] <= 0.0)
240 0 : _tri_sch_mesh._gij_map[iz][i_gap] = 0.0;
241 : }
242 : }
243 5 : }
244 :
245 3331 : for (unsigned int iz = 1; iz < _n_cells + 1; iz++)
246 : {
247 230220 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
248 : {
249 227040 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
250 227040 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
251 227040 : _mdot_soln->set(node_out, (*_mdot_soln)(node_in));
252 : }
253 : }
254 :
255 : // We must do a global assembly to make sure data is parallel consistent before we do things
256 : // like compute L2 norms
257 151 : _aux->solution().close();
258 151 : }
259 :
260 : Real
261 3240330 : TriSubChannel1PhaseProblem::computeAddedHeatPin(unsigned int i_ch, unsigned int iz) const
262 : {
263 : // Compute axial location of nodes.
264 3240330 : auto z2 = _z_grid[iz];
265 3240330 : auto z1 = _z_grid[iz - 1];
266 3240330 : auto heated_length = _subchannel_mesh.getHeatedLength();
267 3240330 : auto unheated_length_entry = _subchannel_mesh.getHeatedLengthEntry();
268 3240330 : if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
269 2740194 : MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
270 : {
271 : // Compute the height of this element.
272 1910652 : auto dz = z2 - z1;
273 1910652 : if (_pin_mesh_exist)
274 : {
275 : double factor;
276 1020936 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
277 1020936 : switch (subch_type)
278 : {
279 : case EChannelType::CENTER:
280 : factor = 1.0 / 6.0;
281 : break;
282 205344 : case EChannelType::EDGE:
283 : factor = 1.0 / 4.0;
284 205344 : break;
285 : case EChannelType::CORNER:
286 : factor = 1.0 / 6.0;
287 : break;
288 : default:
289 : return 0.0; // handle invalid subch_type if needed
290 : }
291 : double heat_rate_in = 0.0;
292 : double heat_rate_out = 0.0;
293 3759408 : for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
294 : {
295 2738472 : auto * node_in = _subchannel_mesh.getPinNode(i_pin, iz - 1);
296 2738472 : auto * node_out = _subchannel_mesh.getPinNode(i_pin, iz);
297 2738472 : heat_rate_out += factor * (*_q_prime_soln)(node_out);
298 2738472 : heat_rate_in += factor * (*_q_prime_soln)(node_in);
299 : }
300 1020936 : return (heat_rate_in + heat_rate_out) * dz / 2.0;
301 : }
302 : else
303 : {
304 889716 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
305 889716 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
306 889716 : return ((*_q_prime_soln)(node_out) + (*_q_prime_soln)(node_in)) * dz / 2.0;
307 : }
308 : }
309 : else
310 : return 0.0;
311 : }
312 :
313 : Real
314 279600 : TriSubChannel1PhaseProblem::getSubChannelPeripheralDuctWidth(unsigned int i_ch) const
315 : {
316 279600 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
317 279600 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
318 : {
319 279600 : auto width = _subchannel_mesh.getPitch();
320 279600 : if (subch_type == EChannelType::CORNER)
321 57600 : width = 2.0 / std::sqrt(3.0) *
322 57600 : (_subchannel_mesh.getPinDiameter() / 2.0 + _tri_sch_mesh.getDuctToPinGap());
323 279600 : return width;
324 : }
325 : else
326 0 : mooseError("Channel is not a perimetric subchannel ");
327 : }
328 :
329 : void
330 1581 : TriSubChannel1PhaseProblem::computeh(int iblock)
331 : {
332 1581 : unsigned int last_node = (iblock + 1) * _block_size;
333 1581 : unsigned int first_node = iblock * _block_size + 1;
334 1581 : const Real & wire_lead_length = _tri_sch_mesh.getWireLeadLength();
335 : const Real & wire_diameter = _tri_sch_mesh.getWireDiameter();
336 1581 : const Real & pitch = _subchannel_mesh.getPitch();
337 1581 : const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
338 :
339 1581 : if (iblock == 0)
340 : {
341 115203 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
342 : {
343 113622 : auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
344 113622 : auto h_out = _fp->h_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
345 113622 : if (h_out < 0)
346 : {
347 0 : mooseError(
348 0 : name(), " : Calculation of negative Enthalpy h_out = : ", h_out, " Axial Level= : ", 0);
349 : }
350 113622 : _h_soln->set(node, h_out);
351 : }
352 : }
353 :
354 1581 : if (!_implicit_bool)
355 : {
356 90 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
357 : {
358 75 : auto z_grid = _subchannel_mesh.getZGrid();
359 75 : auto dz = z_grid[iz] - z_grid[iz - 1];
360 : // Calculation of average mass flux of all periphery subchannels
361 : Real edge_flux_ave = 0.0;
362 : Real mdot_sum = 0.0;
363 : Real si_sum = 0.0;
364 3225 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
365 : {
366 3150 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
367 3150 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
368 : {
369 1350 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
370 1350 : auto Si = (*_S_flow_soln)(node_in);
371 1350 : auto mdot_in = (*_mdot_soln)(node_in);
372 1350 : mdot_sum = mdot_sum + mdot_in;
373 1350 : si_sum = si_sum + Si;
374 : }
375 : }
376 75 : edge_flux_ave = mdot_sum / si_sum;
377 :
378 3225 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
379 : {
380 3150 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
381 3150 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
382 3150 : auto mdot_in = (*_mdot_soln)(node_in);
383 3150 : auto h_in = (*_h_soln)(node_in); // J/kg
384 3150 : auto volume = dz * (*_S_flow_soln)(node_in);
385 3150 : auto mdot_out = (*_mdot_soln)(node_out);
386 3150 : auto h_out = 0.0;
387 : Real sumWijh = 0.0;
388 : Real sumWijPrimeDhij = 0.0;
389 : Real sweep_enthalpy = 0.0;
390 : Real e_cond = 0.0;
391 :
392 : // Calculate added enthalpy from heatflux (Pin, Duct)
393 3150 : Real added_enthalpy = computeAddedHeatPin(i_ch, iz);
394 3150 : added_enthalpy += computeAddedHeatDuct(i_ch, iz);
395 :
396 : // Calculate net sum of enthalpy into/out-of channel i from channels j around i
397 : // (Turbulent diffusion, Diversion Crossflow, Sweep flow Enthalpy, Radial heat conduction)
398 : unsigned int counter = 0;
399 12150 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
400 : {
401 9000 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
402 9000 : auto gap = _subchannel_mesh.getGapWidth(iz, i_gap);
403 9000 : auto Sij = dz * gap;
404 : unsigned int ii_ch = chans.first; // the first subchannel next to gap i_gap
405 : unsigned int jj_ch = chans.second; // the second subchannel next to gap i_gap
406 9000 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
407 9000 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
408 9000 : auto subch_type_i = _subchannel_mesh.getSubchannelType(ii_ch);
409 9000 : auto subch_type_j = _subchannel_mesh.getSubchannelType(jj_ch);
410 : // Define donor enthalpy
411 : auto h_star = 0.0;
412 9000 : if (_Wij(i_gap, iz) > 0.0)
413 7650 : h_star = (*_h_soln)(node_in_i);
414 1350 : else if (_Wij(i_gap, iz) < 0.0)
415 1350 : h_star = (*_h_soln)(node_in_j);
416 : // Diversion crossflow
417 : // take care of the sign by applying the map, use donor cell
418 9000 : sumWijh += _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * h_star;
419 9000 : counter++;
420 : // SWEEP FLOW is calculated if i_gap is located in the periphery
421 : // and we have a wire-wrap (if i_gap is in the periphery then i_ch is in the periphery)
422 : // There are two gaps per periphery subchannel that this is true.
423 9000 : if ((subch_type_i == EChannelType::CORNER || subch_type_i == EChannelType::EDGE) &&
424 2700 : (subch_type_j == EChannelType::CORNER || subch_type_j == EChannelType::EDGE) &&
425 2700 : (wire_lead_length != 0) && (wire_diameter != 0))
426 : {
427 : // donor subchannel and node of sweep flow. The donor subchannel is the subchannel next
428 : // to i_ch that sweep flow, flows from and into i_ch
429 2700 : auto sweep_donor = _tri_sch_mesh.getSweepFlowChans(i_ch).first;
430 2700 : auto * node_sweep_donor = _subchannel_mesh.getChannelNode(sweep_donor, iz - 1);
431 : // if one of the neighbor subchannels of the periphery gap is the donor subchannel
432 : //(the other would be the i_ch) sweep enthalpy flows into i_ch
433 2700 : if ((ii_ch == sweep_donor) || (jj_ch == sweep_donor))
434 : {
435 1350 : sweep_enthalpy += computeSweepFlowMixingParameter(i_gap, iz) * edge_flux_ave * Sij *
436 1350 : (*_h_soln)(node_sweep_donor);
437 : }
438 : // else sweep enthalpy flows out of i_ch
439 : else
440 : {
441 1350 : sweep_enthalpy -= computeSweepFlowMixingParameter(i_gap, iz) * edge_flux_ave * Sij *
442 1350 : (*_h_soln)(node_in);
443 : }
444 : }
445 : // Inner gap
446 : // Turbulent Diffusion
447 : else
448 : {
449 6300 : sumWijPrimeDhij +=
450 6300 : _WijPrime(i_gap, iz) * (2 * h_in - (*_h_soln)(node_in_j) - (*_h_soln)(node_in_i));
451 : }
452 :
453 : // compute the radial heat conduction through the gaps
454 : Real dist_ij = pitch;
455 :
456 9000 : if (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::EDGE)
457 : {
458 900 : dist_ij = pitch;
459 : }
460 8100 : else if ((subch_type_i == EChannelType::CORNER && subch_type_j == EChannelType::EDGE) ||
461 7950 : (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::CORNER))
462 : {
463 1800 : dist_ij = pitch;
464 : }
465 : else
466 : {
467 6300 : dist_ij = pitch / std::sqrt(3);
468 : }
469 :
470 9000 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
471 9000 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
472 : auto shape_factor =
473 9000 : 0.66 * (pitch / pin_diameter) *
474 9000 : std::pow((_subchannel_mesh.getGapWidth(iz, i_gap) / pin_diameter), -0.3);
475 9000 : if (ii_ch == i_ch)
476 : {
477 4500 : e_cond += 0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
478 4500 : ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
479 : }
480 : else
481 : {
482 4500 : e_cond += -0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
483 4500 : ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
484 : }
485 : }
486 :
487 : // compute the axial heat conduction between current and lower axial node
488 3150 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz);
489 3150 : auto * node_in_j = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
490 3150 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
491 3150 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
492 3150 : auto Si = (*_S_flow_soln)(node_in_i);
493 3150 : auto dist_ij = z_grid[iz] - z_grid[iz - 1];
494 :
495 3150 : e_cond += 0.5 * (thcon_i + thcon_j) * Si * ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) /
496 : dist_ij;
497 :
498 3150 : unsigned int nz = _subchannel_mesh.getNumOfAxialCells();
499 : // compute the axial heat conduction between current and upper axial node
500 3150 : if (iz < nz)
501 : {
502 2520 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz);
503 2520 : auto * node_in_j = _subchannel_mesh.getChannelNode(i_ch, iz + 1);
504 2520 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
505 2520 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
506 2520 : auto Si = (*_S_flow_soln)(node_in_i);
507 2520 : auto dist_ij = z_grid[iz + 1] - z_grid[iz];
508 2520 : e_cond += 0.5 * (thcon_i + thcon_j) * Si *
509 2520 : ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
510 : }
511 :
512 : // end of radial heat conduction calc.
513 3150 : h_out =
514 6300 : (mdot_in * h_in - sumWijh - sumWijPrimeDhij + added_enthalpy + e_cond + sweep_enthalpy +
515 3150 : _TR * _rho_soln->old(node_out) * _h_soln->old(node_out) * volume / _dt) /
516 3150 : (mdot_out + _TR * (*_rho_soln)(node_out)*volume / _dt);
517 3150 : if (h_out < 0)
518 : {
519 0 : mooseError(name(),
520 : " : Calculation of negative Enthalpy h_out = : ",
521 : h_out,
522 : " Axial Level= : ",
523 : iz);
524 : }
525 3150 : _h_soln->set(node_out, h_out); // J/kg
526 : }
527 75 : }
528 : }
529 : else
530 : {
531 1566 : LibmeshPetscCall(MatZeroEntries(_hc_time_derivative_mat));
532 1566 : LibmeshPetscCall(MatZeroEntries(_hc_advective_derivative_mat));
533 1566 : LibmeshPetscCall(MatZeroEntries(_hc_cross_derivative_mat));
534 1566 : LibmeshPetscCall(MatZeroEntries(_hc_axial_heat_conduction_mat));
535 1566 : LibmeshPetscCall(MatZeroEntries(_hc_radial_heat_conduction_mat));
536 1566 : LibmeshPetscCall(MatZeroEntries(_hc_sweep_enthalpy_mat));
537 :
538 1566 : LibmeshPetscCall(VecZeroEntries(_hc_time_derivative_rhs));
539 1566 : LibmeshPetscCall(VecZeroEntries(_hc_advective_derivative_rhs));
540 1566 : LibmeshPetscCall(VecZeroEntries(_hc_cross_derivative_rhs));
541 1566 : LibmeshPetscCall(VecZeroEntries(_hc_added_heat_rhs));
542 1566 : LibmeshPetscCall(VecZeroEntries(_hc_axial_heat_conduction_rhs));
543 1566 : LibmeshPetscCall(VecZeroEntries(_hc_radial_heat_conduction_rhs));
544 1566 : LibmeshPetscCall(VecZeroEntries(_hc_sweep_enthalpy_rhs));
545 :
546 1566 : LibmeshPetscCall(MatZeroEntries(_hc_sys_h_mat));
547 1566 : LibmeshPetscCall(VecZeroEntries(_hc_sys_h_rhs));
548 :
549 51606 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
550 : {
551 50040 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
552 50040 : auto pitch = _subchannel_mesh.getPitch();
553 50040 : auto pin_diameter = _subchannel_mesh.getPinDiameter();
554 50040 : auto iz_ind = iz - first_node;
555 :
556 3257400 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
557 : {
558 3207360 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
559 3207360 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
560 3207360 : auto S_in = (*_S_flow_soln)(node_in);
561 3207360 : auto S_out = (*_S_flow_soln)(node_out);
562 3207360 : auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
563 3207360 : auto volume = dz * S_interp;
564 :
565 : PetscScalar Pe = 0.5;
566 3207360 : if (_interpolation_scheme == 3)
567 : {
568 : // Compute the Peclet number
569 785820 : auto w_perim_in = (*_w_perim_soln)(node_in);
570 785820 : auto w_perim_out = (*_w_perim_soln)(node_out);
571 785820 : auto w_perim_interp = this->computeInterpolatedValue(w_perim_out, w_perim_in, 0.5);
572 785820 : auto K_in = _fp->k_from_p_T((*_P_soln)(node_in) + _P_out, (*_T_soln)(node_in));
573 785820 : auto K_out = _fp->k_from_p_T((*_P_soln)(node_out) + _P_out, (*_T_soln)(node_out));
574 785820 : auto K = this->computeInterpolatedValue(K_out, K_in, 0.5);
575 785820 : auto cp_in = _fp->cp_from_p_T((*_P_soln)(node_in) + _P_out, (*_T_soln)(node_in));
576 785820 : auto cp_out = _fp->cp_from_p_T((*_P_soln)(node_out) + _P_out, (*_T_soln)(node_out));
577 785820 : auto cp = this->computeInterpolatedValue(cp_out, cp_in, 0.5);
578 : auto mdot_loc =
579 785820 : this->computeInterpolatedValue((*_mdot_soln)(node_out), (*_mdot_soln)(node_in), 0.5);
580 : // hydraulic diameter in the i direction
581 785820 : auto Dh_i = 4.0 * S_interp / w_perim_interp;
582 785820 : Pe = mdot_loc * Dh_i * cp / (K * S_interp) * (mdot_loc / std::abs(mdot_loc));
583 : }
584 3207360 : auto alpha = computeInterpolationCoefficients(Pe);
585 :
586 : // Time derivative term
587 3207360 : if (iz == first_node)
588 : {
589 : PetscScalar value_vec_tt =
590 112992 : -1.0 * _TR * alpha * (*_rho_soln)(node_in) * (*_h_soln)(node_in)*volume / _dt;
591 112992 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
592 112992 : LibmeshPetscCall(
593 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
594 : }
595 : else
596 : {
597 3094368 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
598 3094368 : PetscInt col_tt = i_ch + _n_channels * (iz_ind - 1);
599 3094368 : PetscScalar value_tt = _TR * alpha * (*_rho_soln)(node_in)*volume / _dt;
600 3094368 : LibmeshPetscCall(MatSetValues(
601 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
602 : }
603 :
604 : // Adding diagonal elements
605 3207360 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
606 3207360 : PetscInt col_tt = i_ch + _n_channels * iz_ind;
607 3207360 : PetscScalar value_tt = _TR * (1.0 - alpha) * (*_rho_soln)(node_out)*volume / _dt;
608 3207360 : LibmeshPetscCall(MatSetValues(
609 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
610 :
611 : // Adding RHS elements
612 : PetscScalar rho_old_interp =
613 3207360 : computeInterpolatedValue(_rho_soln->old(node_out), _rho_soln->old(node_in), Pe);
614 : PetscScalar h_old_interp =
615 3207360 : computeInterpolatedValue(_h_soln->old(node_out), _h_soln->old(node_in), Pe);
616 3207360 : PetscScalar value_vec_tt = _TR * rho_old_interp * h_old_interp * volume / _dt;
617 3207360 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
618 3207360 : LibmeshPetscCall(
619 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
620 :
621 : // Advective derivative term
622 3207360 : if (iz == first_node)
623 : {
624 112992 : PetscInt row_at = i_ch + _n_channels * iz_ind;
625 112992 : PetscScalar value_at = alpha * (*_mdot_soln)(node_in) * (*_h_soln)(node_in);
626 112992 : LibmeshPetscCall(
627 : VecSetValues(_hc_advective_derivative_rhs, 1, &row_at, &value_at, ADD_VALUES));
628 :
629 112992 : value_at = alpha * (*_mdot_soln)(node_out) - (1 - alpha) * (*_mdot_soln)(node_in);
630 112992 : PetscInt col_at = i_ch + _n_channels * iz_ind;
631 112992 : LibmeshPetscCall(MatSetValues(
632 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
633 :
634 112992 : value_at = (1 - alpha) * (*_mdot_soln)(node_out);
635 112992 : col_at = i_ch + _n_channels * (iz_ind + 1);
636 112992 : LibmeshPetscCall(MatSetValues(
637 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
638 : }
639 3094368 : else if (iz == last_node)
640 : {
641 112992 : PetscInt row_at = i_ch + _n_channels * iz_ind;
642 112992 : PetscScalar value_at = 1.0 * (*_mdot_soln)(node_out);
643 112992 : PetscInt col_at = i_ch + _n_channels * iz_ind;
644 112992 : LibmeshPetscCall(MatSetValues(
645 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
646 :
647 112992 : value_at = -1.0 * (*_mdot_soln)(node_in);
648 112992 : col_at = i_ch + _n_channels * (iz_ind - 1);
649 112992 : LibmeshPetscCall(MatSetValues(
650 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
651 : }
652 : else
653 : {
654 2981376 : PetscInt row_at = i_ch + _n_channels * iz_ind;
655 : PetscInt col_at;
656 :
657 2981376 : PetscScalar value_at = -alpha * (*_mdot_soln)(node_in);
658 2981376 : col_at = i_ch + _n_channels * (iz_ind - 1);
659 2981376 : LibmeshPetscCall(MatSetValues(
660 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
661 :
662 2981376 : value_at = alpha * (*_mdot_soln)(node_out) - (1 - alpha) * (*_mdot_soln)(node_in);
663 2981376 : col_at = i_ch + _n_channels * iz_ind;
664 2981376 : LibmeshPetscCall(MatSetValues(
665 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
666 :
667 2981376 : value_at = (1 - alpha) * (*_mdot_soln)(node_out);
668 2981376 : col_at = i_ch + _n_channels * (iz_ind + 1);
669 2981376 : LibmeshPetscCall(MatSetValues(
670 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, ADD_VALUES));
671 : }
672 :
673 : // Axial heat conduction
674 3207360 : auto * node_center = _subchannel_mesh.getChannelNode(i_ch, iz);
675 3207360 : auto K_center = _fp->k_from_p_T((*_P_soln)(node_center) + _P_out, (*_T_soln)(node_center));
676 : auto cp_center =
677 3207360 : _fp->cp_from_p_T((*_P_soln)(node_center) + _P_out, (*_T_soln)(node_center));
678 3207360 : auto diff_center = K_center / (cp_center + 1e-15);
679 :
680 3207360 : if (iz == first_node)
681 : {
682 112992 : auto * node_top = _subchannel_mesh.getChannelNode(i_ch, iz + 1);
683 112992 : auto * node_bottom = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
684 : auto K_bottom =
685 112992 : _fp->k_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
686 112992 : auto K_top = _fp->k_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
687 : auto cp_bottom =
688 112992 : _fp->cp_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
689 112992 : auto cp_top = _fp->cp_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
690 112992 : auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
691 112992 : auto diff_top = K_top / (cp_top + 1e-15);
692 :
693 112992 : auto dz_up = _z_grid[iz + 1] - _z_grid[iz];
694 112992 : auto dz_down = _z_grid[iz] - _z_grid[iz - 1];
695 : auto S_up =
696 112992 : computeInterpolatedValue((*_S_flow_soln)(node_top), (*_S_flow_soln)(node_center));
697 : auto S_down =
698 112992 : computeInterpolatedValue((*_S_flow_soln)(node_center), (*_S_flow_soln)(node_bottom));
699 112992 : auto diff_up = computeInterpolatedValue(diff_top, diff_center);
700 112992 : auto diff_down = computeInterpolatedValue(diff_center, diff_bottom);
701 :
702 : // Diagonal value
703 112992 : PetscInt row_at = i_ch + _n_channels * iz_ind;
704 112992 : PetscInt col_at = i_ch + _n_channels * iz_ind;
705 112992 : PetscScalar value_at = diff_up * S_up / dz_up + diff_down * S_down / dz_down;
706 112992 : LibmeshPetscCall(MatSetValues(
707 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
708 :
709 : // Bottom value
710 112992 : value_at = 1.0 * diff_down * S_down / dz_down * (*_h_soln)(node_bottom);
711 112992 : LibmeshPetscCall(
712 : VecSetValues(_hc_axial_heat_conduction_rhs, 1, &row_at, &value_at, ADD_VALUES));
713 :
714 : // Top value
715 112992 : col_at = i_ch + _n_channels * (iz_ind + 1);
716 112992 : value_at = -diff_up * S_up / dz_up;
717 112992 : LibmeshPetscCall(MatSetValues(
718 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
719 : }
720 3094368 : else if (iz == last_node)
721 : {
722 112992 : auto * node_bottom = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
723 : auto K_bottom =
724 112992 : _fp->k_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
725 : auto cp_bottom =
726 112992 : _fp->cp_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
727 112992 : auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
728 :
729 112992 : auto dz_down = _z_grid[iz] - _z_grid[iz - 1];
730 112992 : auto S_down = 0.5 * ((*_S_flow_soln)(node_center) + (*_S_flow_soln)(node_bottom));
731 112992 : auto diff_down = 0.5 * (diff_center + diff_bottom);
732 :
733 : // Diagonal value
734 112992 : PetscInt row_at = i_ch + _n_channels * iz_ind;
735 112992 : PetscInt col_at = i_ch + _n_channels * iz_ind;
736 112992 : PetscScalar value_at = diff_down * S_down / dz_down;
737 112992 : LibmeshPetscCall(MatSetValues(
738 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
739 :
740 : // Bottom value
741 112992 : col_at = i_ch + _n_channels * (iz_ind - 1);
742 112992 : value_at = -diff_down * S_down / dz_down;
743 112992 : LibmeshPetscCall(MatSetValues(
744 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
745 :
746 : // Outflow derivative
747 : /// TODO: Current axial derivative is zero - check if outflow conditions may make a difference
748 : // value_at = -1.0 * (*_mdot_soln)(node_center) * (*_h_soln)(node_center);
749 : // VecSetValues(_hc_axial_heat_conduction_rhs, 1, &row_at, &value_at, ADD_VALUES);
750 : }
751 : else
752 : {
753 2981376 : auto * node_top = _subchannel_mesh.getChannelNode(i_ch, iz + 1);
754 2981376 : auto * node_bottom = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
755 : auto K_bottom =
756 2981376 : _fp->k_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
757 2981376 : auto K_top = _fp->k_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
758 : auto cp_bottom =
759 2981376 : _fp->cp_from_p_T((*_P_soln)(node_bottom) + _P_out, (*_T_soln)(node_bottom));
760 2981376 : auto cp_top = _fp->cp_from_p_T((*_P_soln)(node_top) + _P_out, (*_T_soln)(node_top));
761 2981376 : auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
762 2981376 : auto diff_top = K_top / (cp_top + 1e-15);
763 :
764 2981376 : auto dz_up = _z_grid[iz + 1] - _z_grid[iz];
765 2981376 : auto dz_down = _z_grid[iz] - _z_grid[iz - 1];
766 : auto S_up =
767 2981376 : computeInterpolatedValue((*_S_flow_soln)(node_top), (*_S_flow_soln)(node_center));
768 : auto S_down =
769 2981376 : computeInterpolatedValue((*_S_flow_soln)(node_center), (*_S_flow_soln)(node_bottom));
770 2981376 : auto diff_up = computeInterpolatedValue(diff_top, diff_center);
771 2981376 : auto diff_down = computeInterpolatedValue(diff_center, diff_bottom);
772 :
773 : // Diagonal value
774 2981376 : PetscInt row_at = i_ch + _n_channels * iz_ind;
775 2981376 : PetscInt col_at = i_ch + _n_channels * iz_ind;
776 2981376 : PetscScalar value_at = diff_up * S_up / dz_up + diff_down * S_down / dz_down;
777 2981376 : LibmeshPetscCall(MatSetValues(
778 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
779 :
780 : // Bottom value
781 2981376 : col_at = i_ch + _n_channels * (iz_ind - 1);
782 2981376 : value_at = -diff_down * S_down / dz_down;
783 2981376 : LibmeshPetscCall(MatSetValues(
784 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
785 :
786 : // Top value
787 2981376 : col_at = i_ch + _n_channels * (iz_ind + 1);
788 2981376 : value_at = -diff_up * S_up / dz_up;
789 2981376 : LibmeshPetscCall(MatSetValues(
790 : _hc_axial_heat_conduction_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
791 : }
792 :
793 : // Radial Terms
794 : unsigned int counter = 0;
795 : unsigned int cross_index = iz;
796 : // Real radial_heat_conduction(0.0);
797 12529200 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
798 : {
799 9321840 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
800 : unsigned int ii_ch = chans.first;
801 : unsigned int jj_ch = chans.second;
802 9321840 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
803 9321840 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
804 : PetscScalar h_star;
805 : // figure out donor axial velocity
806 9321840 : if (_Wij(i_gap, cross_index) > 0.0)
807 : {
808 6566710 : if (iz == first_node)
809 : {
810 245094 : h_star = (*_h_soln)(node_in_i);
811 490188 : PetscScalar value_vec_ct = -1.0 * alpha *
812 245094 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
813 245094 : _Wij(i_gap, cross_index) * h_star;
814 245094 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
815 245094 : LibmeshPetscCall(VecSetValues(
816 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
817 : }
818 : else
819 : {
820 6321616 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
821 6321616 : _Wij(i_gap, cross_index);
822 6321616 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
823 6321616 : PetscInt col_ct = ii_ch + _n_channels * (iz_ind - 1);
824 6321616 : LibmeshPetscCall(MatSetValues(
825 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
826 : }
827 6566710 : PetscScalar value_ct = (1.0 - alpha) *
828 6566710 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
829 6566710 : _Wij(i_gap, cross_index);
830 6566710 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
831 6566710 : PetscInt col_ct = ii_ch + _n_channels * iz_ind;
832 6566710 : LibmeshPetscCall(MatSetValues(
833 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
834 : }
835 2755130 : else if (_Wij(i_gap, cross_index) < 0.0) // _Wij=0 operations not necessary
836 : {
837 2677658 : if (iz == first_node)
838 : {
839 80754 : h_star = (*_h_soln)(node_in_j);
840 161508 : PetscScalar value_vec_ct = -1.0 * alpha *
841 80754 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
842 80754 : _Wij(i_gap, cross_index) * h_star;
843 80754 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
844 80754 : LibmeshPetscCall(VecSetValues(
845 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
846 : }
847 : else
848 : {
849 2596904 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
850 2596904 : _Wij(i_gap, cross_index);
851 2596904 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
852 2596904 : PetscInt col_ct = jj_ch + _n_channels * (iz_ind - 1);
853 2596904 : LibmeshPetscCall(MatSetValues(
854 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
855 : }
856 2677658 : PetscScalar value_ct = (1.0 - alpha) *
857 2677658 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
858 2677658 : _Wij(i_gap, cross_index);
859 2677658 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
860 2677658 : PetscInt col_ct = jj_ch + _n_channels * iz_ind;
861 2677658 : LibmeshPetscCall(MatSetValues(
862 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
863 : }
864 :
865 : // Turbulent cross flows
866 9321840 : if (iz == first_node)
867 : {
868 : PetscScalar value_vec_ct =
869 329580 : -2.0 * alpha * (*_h_soln)(node_in)*_WijPrime(i_gap, cross_index);
870 329580 : value_vec_ct += alpha * (*_h_soln)(node_in_j)*_WijPrime(i_gap, cross_index);
871 329580 : value_vec_ct += alpha * (*_h_soln)(node_in_i)*_WijPrime(i_gap, cross_index);
872 329580 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
873 329580 : LibmeshPetscCall(
874 : VecSetValues(_hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
875 : }
876 : else
877 : {
878 8992260 : PetscScalar value_center_ct = 2.0 * alpha * _WijPrime(i_gap, cross_index);
879 8992260 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
880 8992260 : PetscInt col_ct = i_ch + _n_channels * (iz_ind - 1);
881 8992260 : LibmeshPetscCall(MatSetValues(
882 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
883 :
884 8992260 : PetscScalar value_left_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
885 8992260 : row_ct = i_ch + _n_channels * iz_ind;
886 8992260 : col_ct = jj_ch + _n_channels * (iz_ind - 1);
887 8992260 : LibmeshPetscCall(MatSetValues(
888 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
889 :
890 8992260 : PetscScalar value_right_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
891 8992260 : row_ct = i_ch + _n_channels * iz_ind;
892 8992260 : col_ct = ii_ch + _n_channels * (iz_ind - 1);
893 8992260 : LibmeshPetscCall(MatSetValues(
894 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
895 : }
896 9321840 : PetscScalar value_center_ct = 2.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
897 9321840 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
898 9321840 : PetscInt col_ct = i_ch + _n_channels * iz_ind;
899 9321840 : LibmeshPetscCall(MatSetValues(
900 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
901 :
902 9321840 : PetscScalar value_left_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
903 9321840 : row_ct = i_ch + _n_channels * iz_ind;
904 9321840 : col_ct = jj_ch + _n_channels * iz_ind;
905 9321840 : LibmeshPetscCall(MatSetValues(
906 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
907 :
908 9321840 : PetscScalar value_right_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
909 9321840 : row_ct = i_ch + _n_channels * iz_ind;
910 9321840 : col_ct = ii_ch + _n_channels * iz_ind;
911 9321840 : LibmeshPetscCall(MatSetValues(
912 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
913 :
914 : // Radial heat conduction
915 9321840 : auto subch_type_i = _subchannel_mesh.getSubchannelType(ii_ch);
916 9321840 : auto subch_type_j = _subchannel_mesh.getSubchannelType(jj_ch);
917 : Real dist_ij = pitch;
918 :
919 9321840 : if (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::EDGE)
920 : {
921 : dist_ij = pitch;
922 : }
923 8401200 : else if ((subch_type_i == EChannelType::CORNER && subch_type_j == EChannelType::EDGE) ||
924 8301120 : (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::CORNER))
925 : {
926 : dist_ij = pitch;
927 : }
928 : else
929 : {
930 7200240 : dist_ij = pitch / std::sqrt(3);
931 : }
932 :
933 9321840 : auto Sij = dz * _subchannel_mesh.getGapWidth(iz, i_gap);
934 9321840 : auto K_i = _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
935 9321840 : auto K_j = _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
936 9321840 : auto cp_i = _fp->cp_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i));
937 9321840 : auto cp_j = _fp->cp_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j));
938 9321840 : auto A_i = K_i / cp_i;
939 9321840 : auto A_j = K_j / cp_j;
940 9321840 : auto harm_A = 2.0 * A_i * A_j / (A_i + A_j);
941 : auto shape_factor =
942 9321840 : 0.66 * (pitch / pin_diameter) *
943 9321840 : std::pow((_subchannel_mesh.getGapWidth(iz, i_gap) / pin_diameter), -0.3);
944 : // auto base_value = 0.5 * (A_i + A_j) * Sij * shape_factor / dist_ij;
945 9321840 : auto base_value = harm_A * shape_factor * Sij / dist_ij;
946 9321840 : auto neg_base_value = -1.0 * base_value;
947 :
948 9321840 : row_ct = ii_ch + _n_channels * iz_ind;
949 9321840 : col_ct = ii_ch + _n_channels * iz_ind;
950 9321840 : LibmeshPetscCall(MatSetValues(
951 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &base_value, ADD_VALUES));
952 :
953 9321840 : row_ct = jj_ch + _n_channels * iz_ind;
954 9321840 : col_ct = jj_ch + _n_channels * iz_ind;
955 9321840 : LibmeshPetscCall(MatSetValues(
956 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &base_value, ADD_VALUES));
957 :
958 9321840 : row_ct = ii_ch + _n_channels * iz_ind;
959 9321840 : col_ct = jj_ch + _n_channels * iz_ind;
960 9321840 : LibmeshPetscCall(MatSetValues(
961 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &neg_base_value, ADD_VALUES));
962 :
963 9321840 : row_ct = jj_ch + _n_channels * iz_ind;
964 9321840 : col_ct = ii_ch + _n_channels * iz_ind;
965 9321840 : LibmeshPetscCall(MatSetValues(
966 : _hc_radial_heat_conduction_mat, 1, &row_ct, 1, &col_ct, &neg_base_value, ADD_VALUES));
967 9321840 : counter++;
968 : }
969 :
970 : // Compute the sweep flow enthalpy change
971 : // Calculation of average mass flux of all periphery subchannels
972 : Real edge_flux_ave = 0.0;
973 : Real mdot_sum = 0.0;
974 : Real si_sum = 0.0;
975 272946720 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
976 : {
977 269739360 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
978 269739360 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
979 : {
980 77188320 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
981 77188320 : auto Si = (*_S_flow_soln)(node_in);
982 77188320 : auto mdot_in = (*_mdot_soln)(node_in);
983 77188320 : mdot_sum = mdot_sum + mdot_in;
984 77188320 : si_sum = si_sum + Si;
985 : }
986 : }
987 3207360 : edge_flux_ave = mdot_sum / si_sum;
988 3207360 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
989 : PetscScalar sweep_enthalpy = 0.0;
990 3207360 : if ((subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER) &&
991 1060800 : (wire_diameter != 0.0) && (wire_lead_length != 0.0))
992 : {
993 : auto beta_in = std::numeric_limits<double>::quiet_NaN();
994 : auto beta_out = std::numeric_limits<double>::quiet_NaN();
995 : // donor sweep channel for i_ch
996 970800 : auto sweep_donor = _tri_sch_mesh.getSweepFlowChans(i_ch).first;
997 970800 : auto * node_sweep_donor = _subchannel_mesh.getChannelNode(sweep_donor, iz - 1);
998 : // Calculation of turbulent mixing parameter
999 3612960 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
1000 : {
1001 2642160 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
1002 : unsigned int ii_ch = chans.first;
1003 : unsigned int jj_ch = chans.second;
1004 2642160 : auto subch_type_i = _subchannel_mesh.getSubchannelType(ii_ch);
1005 2642160 : auto subch_type_j = _subchannel_mesh.getSubchannelType(jj_ch);
1006 2642160 : if ((subch_type_i == EChannelType::CORNER || subch_type_i == EChannelType::EDGE) &&
1007 1941600 : (subch_type_j == EChannelType::CORNER || subch_type_j == EChannelType::EDGE))
1008 : {
1009 1941600 : if ((ii_ch == sweep_donor) || (jj_ch == sweep_donor))
1010 : {
1011 970800 : beta_in = computeSweepFlowMixingParameter(i_gap, iz);
1012 : }
1013 : else
1014 : {
1015 970800 : beta_out = computeSweepFlowMixingParameter(i_gap, iz);
1016 : }
1017 : }
1018 : }
1019 : // Abort execution if required values are unset
1020 : mooseAssert(!std::isnan(beta_in),
1021 : "beta_in was not set. Check gap logic for i_ch = " + std::to_string(i_ch) +
1022 : ", iz = " + std::to_string(iz));
1023 : mooseAssert(!std::isnan(beta_out),
1024 : "beta_out was not set. Check gap logic for i_ch = " + std::to_string(i_ch) +
1025 : ", iz = " + std::to_string(iz));
1026 :
1027 970800 : auto gap = _tri_sch_mesh.getDuctToPinGap();
1028 970800 : auto Sij = dz * gap;
1029 970800 : auto wsweep_in = edge_flux_ave * beta_in * Sij;
1030 970800 : auto wsweep_out = edge_flux_ave * beta_out * Sij;
1031 970800 : auto sweep_hin = (*_h_soln)(node_sweep_donor);
1032 970800 : auto sweep_hout = (*_h_soln)(node_in);
1033 970800 : sweep_enthalpy = (wsweep_in * sweep_hin - wsweep_out * sweep_hout);
1034 :
1035 970800 : if (iz == first_node)
1036 : {
1037 33240 : PetscInt row_sh = i_ch + _n_channels * iz_ind;
1038 33240 : PetscScalar value_hs = -sweep_enthalpy;
1039 33240 : LibmeshPetscCall(
1040 : VecSetValues(_hc_sweep_enthalpy_rhs, 1, &row_sh, &value_hs, ADD_VALUES));
1041 : }
1042 : else
1043 : {
1044 : // coefficient of sweep_hin
1045 937560 : PetscInt row_sh = i_ch + _n_channels * (iz_ind - 1);
1046 937560 : PetscInt col_sh = i_ch + _n_channels * (iz_ind - 1);
1047 937560 : LibmeshPetscCall(MatSetValues(
1048 : _hc_sweep_enthalpy_mat, 1, &row_sh, 1, &col_sh, &wsweep_out, ADD_VALUES));
1049 937560 : PetscInt col_sh_l = sweep_donor + _n_channels * (iz_ind - 1);
1050 937560 : PetscScalar neg_sweep_in = -1.0 * wsweep_in;
1051 : // coefficient of sweep_hout
1052 937560 : LibmeshPetscCall(MatSetValues(
1053 : _hc_sweep_enthalpy_mat, 1, &row_sh, 1, &col_sh_l, &(neg_sweep_in), ADD_VALUES));
1054 : }
1055 : }
1056 :
1057 : // Add heat enthalpy from pin and/or duct
1058 3207360 : PetscScalar added_enthalpy = computeAddedHeatPin(i_ch, iz);
1059 3207360 : added_enthalpy += computeAddedHeatDuct(i_ch, iz);
1060 3207360 : PetscInt row_vec_ht = i_ch + _n_channels * iz_ind;
1061 3207360 : LibmeshPetscCall(
1062 : VecSetValues(_hc_added_heat_rhs, 1, &row_vec_ht, &added_enthalpy, ADD_VALUES));
1063 : }
1064 : }
1065 : // Assembling system
1066 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
1067 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
1068 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
1069 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
1070 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
1071 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
1072 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_axial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1073 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_axial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1074 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_radial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1075 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_radial_heat_conduction_mat, MAT_FINAL_ASSEMBLY));
1076 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_sweep_enthalpy_mat, MAT_FINAL_ASSEMBLY));
1077 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_sweep_enthalpy_mat, MAT_FINAL_ASSEMBLY));
1078 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1079 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1080 : // Add all matrices together
1081 1566 : LibmeshPetscCall(
1082 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1083 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1084 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1085 1566 : LibmeshPetscCall(
1086 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1087 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1088 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1089 1566 : LibmeshPetscCall(
1090 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1091 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1092 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1093 1566 : LibmeshPetscCall(
1094 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_axial_heat_conduction_mat, DIFFERENT_NONZERO_PATTERN));
1095 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1096 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1097 1566 : LibmeshPetscCall(
1098 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_radial_heat_conduction_mat, DIFFERENT_NONZERO_PATTERN));
1099 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1100 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1101 1566 : LibmeshPetscCall(
1102 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_sweep_enthalpy_mat, DIFFERENT_NONZERO_PATTERN));
1103 1566 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1104 1566 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1105 1566 : if (_verbose_subchannel)
1106 1337 : _console << "Block: " << iblock << " - Enthalpy conservation matrix assembled" << std::endl;
1107 : // RHS
1108 1566 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_time_derivative_rhs));
1109 1566 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_advective_derivative_rhs));
1110 1566 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_cross_derivative_rhs));
1111 1566 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_added_heat_rhs));
1112 1566 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_axial_heat_conduction_rhs));
1113 1566 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_radial_heat_conduction_rhs));
1114 1566 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_sweep_enthalpy_rhs));
1115 :
1116 : // Use system to solve for and populate enthalpy
1117 1566 : LibmeshPetscCall(this->solveAndPopulateEnthalpy(
1118 : _hc_sys_h_mat, _hc_sys_h_rhs, first_node, last_node, "h_sys_"));
1119 : }
1120 1581 : }
|
