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