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