https://mooseframework.inl.gov
QuadSubChannel1PhaseProblem.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 "SCM.h"
13 
15 
18 {
20  params.addClassDescription(
21  "Solver class for subchannels in a square lattice assembly and bare fuel pins");
22  params.addRequiredParam<Real>("beta",
23  "Thermal diffusion coefficient used in turbulent crossflow.");
24  params.addParam<bool>(
25  "default_friction_model",
26  true,
27  "Boolean to define which friction model to use (default: Pang, B. et al. "
28  "KIT, 2013. / non-default: Todreas-Kazimi NUCLEAR SYSTEMS, second edition, Volume 1, 2011)");
29  params.addParam<bool>(
30  "constant_beta",
31  true,
32  "Boolean to define the use of a constant beta or beta correlation (Kim and Chung, 2001)");
33  return params;
34 }
35 
37  : SubChannel1PhaseProblem(params),
38  _subchannel_mesh(SCM::getMesh<QuadSubChannelMesh>(_mesh)),
39  _beta(getParam<Real>("beta")),
40  _default_friction_model(getParam<bool>("default_friction_model")),
41  _constant_beta(getParam<bool>("constant_beta"))
42 {
43 }
44 
45 void
47 {
49  if (_deformation)
50  {
51  Real standard_area, additional_area, wetted_perimeter, displaced_area;
53  auto gap = _subchannel_mesh.getGap();
54  auto z_blockage = _subchannel_mesh.getZBlockage();
55  auto index_blockage = _subchannel_mesh.getIndexBlockage();
56  auto reduction_blockage = _subchannel_mesh.getReductionBlockage();
57  for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
58  {
59  for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
60  {
61  auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
62  auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
63  auto Z = _z_grid[iz];
64  Real rod_area = 0.0;
65  Real rod_perimeter = 0.0;
66  for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
67  {
68  auto * pin_node = _subchannel_mesh.getPinNode(i_pin, iz);
69  rod_area += 0.25 * 0.25 * M_PI * (*_Dpin_soln)(pin_node) * (*_Dpin_soln)(pin_node);
70  rod_perimeter += 0.25 * M_PI * (*_Dpin_soln)(pin_node);
71  }
72 
73  if (subch_type == EChannelType::CORNER)
74  {
75  standard_area = 0.25 * pitch * pitch;
76  displaced_area = (2 * gap + pitch) * (*_displacement_soln)(node) / sqrt(2) +
77  (*_displacement_soln)(node) * (*_displacement_soln)(node) / 2;
78  additional_area = pitch * gap + gap * gap;
79  wetted_perimeter =
80  rod_perimeter + pitch + 2 * gap + 2 * (*_displacement_soln)(node) / sqrt(2);
81  }
82  else if (subch_type == EChannelType::EDGE)
83  {
84  standard_area = 0.5 * pitch * pitch;
85  additional_area = pitch * gap;
86  displaced_area = pitch * (*_displacement_soln)(node);
87  wetted_perimeter = rod_perimeter + pitch;
88  }
89  else
90  {
91  standard_area = pitch * pitch;
92  displaced_area = 0.0;
93  additional_area = 0.0;
94  wetted_perimeter = rod_perimeter;
95  }
96 
98  auto subchannel_area = displaced_area + standard_area + additional_area - rod_area;
99 
101  auto overlapping_pin_area = 0.0;
102  auto overlapping_wetted_perimeter = 0.0;
103  for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
104  {
105  auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
106  auto pin_1 = gap_pins.first;
107  auto pin_2 = gap_pins.second;
108  auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
109  auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
110  auto Diameter1 = (*_Dpin_soln)(pin_node_1);
111  auto Radius1 = Diameter1 / 2.0;
112  auto Diameter2 = (*_Dpin_soln)(pin_node_2);
113  auto Radius2 = Diameter2 / 2.0;
115 
116  if (pitch < (Radius1 + Radius2)) // overlapping pins
117  {
118  mooseWarning(" The gap of index : '", i_gap, " at axial cell ", iz, " ' is blocked.");
119  auto cos1 =
120  (pitch * pitch + Radius1 * Radius1 - Radius2 * Radius2) / (2 * pitch * Radius1);
121  auto cos2 =
122  (pitch * pitch + Radius2 * Radius2 - Radius1 * Radius1) / (2 * pitch * Radius2);
123  auto angle1 = 2.0 * acos(cos1);
124  auto angle2 = 2.0 * acos(cos2);
125  // half of the intersecting arc-length
126  overlapping_wetted_perimeter += 0.5 * angle1 * Radius1 + 0.5 * angle2 * Radius2;
127  // Half of the overlapping area
128  overlapping_pin_area +=
129  0.5 * Radius1 * Radius1 * acos(cos1) + 0.5 * Radius2 * Radius2 * acos(cos2) -
130  0.25 * sqrt((-pitch + Radius1 + Radius2) * (pitch + Radius1 - Radius2) *
131  (pitch - Radius1 + Radius2) * (pitch + Radius1 + Radius2));
132  }
133  }
134  subchannel_area += overlapping_pin_area; // correct surface area
135  wetted_perimeter += -overlapping_wetted_perimeter; // correct wetted perimeter
136 
138  auto index = 0;
139  for (const auto & i_blockage : index_blockage)
140  {
141  if (i_ch == i_blockage && (Z >= z_blockage.front() && Z <= z_blockage.back()))
142  {
143  subchannel_area *= reduction_blockage[index];
144  }
145  index++;
146  }
147 
148  _S_flow_soln->set(node, subchannel_area);
149  _w_perim_soln->set(node, wetted_perimeter);
150  }
151  }
153  for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
154  {
155  for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
156  {
157  auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
158  auto pin_1 = gap_pins.first;
159  auto pin_2 = gap_pins.second;
160  auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
161  auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
162  if (pin_1 == pin_2) // Corner or edge gap
163  {
164  auto displacement = 0.0;
165  auto counter = 0.0;
166  for (auto i_ch : _subchannel_mesh.getPinChannels(pin_1))
167  {
168  auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
169  auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
170  if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
171  {
172  displacement += (*_displacement_soln)(node);
173  counter += 1.0;
174  }
175  }
176  displacement = displacement / counter;
177  _subchannel_mesh._gij_map[iz][i_gap] =
178  (pitch - (*_Dpin_soln)(pin_node_1)) / 2.0 + gap + displacement;
179  }
180  else // center gap
181  {
182  _subchannel_mesh._gij_map[iz][i_gap] =
183  pitch - (*_Dpin_soln)(pin_node_1) / 2.0 - (*_Dpin_soln)(pin_node_2) / 2.0;
184  }
185  // if pins come in contact, the gap is zero
186  if (_subchannel_mesh._gij_map[iz][i_gap] <= 0.0)
187  _subchannel_mesh._gij_map[iz][i_gap] = 0.0;
188  }
189  }
190  }
191 
192  for (unsigned int iz = 1; iz < _n_cells + 1; iz++)
193  {
194  for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
195  {
196  auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
197  auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
198  _mdot_soln->set(node_out, (*_mdot_soln)(node_in));
199  }
200  }
201 
202  // We must do a global assembly to make sure data is parallel consistent before we do things
203  // like compute L2 norms
204  _aux->solution().close();
205 }
206 
207 Real
209 {
210  auto Re = friction_args.Re;
211  auto i_ch = friction_args.i_ch;
214  {
215  Real a, b;
216  if (Re < 1)
217  {
218  return 64.0;
219  }
220  else if (Re >= 1 and Re < 5000)
221  {
222  a = 64.0;
223  b = -1.0;
224  }
225  else if (Re >= 5000 and Re < 30000)
226  {
227  a = 0.316;
228  b = -0.25;
229  }
230  else
231  {
232  a = 0.184;
233  b = -0.20;
234  }
235  return a * std::pow(Re, b);
236  }
238  else
239  {
240  Real aL, b1L, b2L, cL;
241  Real aT, b1T, b2T, cT;
243  auto pin_diameter = _subchannel_mesh.getPinDiameter();
244  // This gap is a constant value for the whole assembly. Might want to make it
245  // subchannel specific in the future if we have duct deformation.
246  auto gap = _subchannel_mesh.getGap();
247  auto w = (pin_diameter / 2.0) + (pitch / 2.0) + gap;
248  auto p_over_d = pitch / pin_diameter;
249  auto w_over_d = w / pin_diameter;
250  auto ReL = std::pow(10, (p_over_d - 1)) * 320.0;
251  auto ReT = std::pow(10, 0.7 * (p_over_d - 1)) * 1.0E+4;
252  auto psi = std::log(Re / ReL) / std::log(ReT / ReL);
253  auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
254  const Real lambda = 7.0;
255 
256  // Find the coefficients of bare Pin bundle friction factor
257  // correlations for turbulent and laminar flow regimes. Todreas & Kazimi, Nuclear Systems Volume
258  // 1
259  if (subch_type == EChannelType::CENTER)
260  {
261  if (p_over_d < 1.1)
262  {
263  aL = 26.37;
264  b1L = 374.2;
265  b2L = -493.9;
266  aT = 0.09423;
267  b1T = 0.5806;
268  b2T = -1.239;
269  }
270  else
271  {
272  aL = 35.55;
273  b1L = 263.7;
274  b2L = -190.2;
275  aT = 0.1339;
276  b1T = 0.09059;
277  b2T = -0.09926;
278  }
279  // laminar flow friction factor for bare Pin bundle - Center subchannel
280  cL = aL + b1L * (p_over_d - 1) + b2L * std::pow((p_over_d - 1), 2);
281  // turbulent flow friction factor for bare Pin bundle - Center subchannel
282  cT = aT + b1T * (p_over_d - 1) + b2T * std::pow((p_over_d - 1), 2);
283  }
284  else if (subch_type == EChannelType::EDGE)
285  {
286  if (p_over_d < 1.1)
287  {
288  aL = 26.18;
289  b1L = 554.5;
290  b2L = -1480;
291  aT = 0.09377;
292  b1T = 0.8732;
293  b2T = -3.341;
294  }
295  else
296  {
297  aL = 44.40;
298  b1L = 256.7;
299  b2L = -267.6;
300  aT = 0.1430;
301  b1T = 0.04199;
302  b2T = -0.04428;
303  }
304  // laminar flow friction factor for bare Pin bundle - Edge subchannel
305  cL = aL + b1L * (w_over_d - 1) + b2L * std::pow((w_over_d - 1), 2);
306  // turbulent flow friction factor for bare Pin bundle - Edge subchannel
307  cT = aT + b1T * (w_over_d - 1) + b2T * std::pow((w_over_d - 1), 2);
308  }
309  else
310  {
311  if (p_over_d < 1.1)
312  {
313  aL = 28.62;
314  b1L = 715.9;
315  b2L = -2807;
316  aT = 0.09755;
317  b1T = 1.127;
318  b2T = -6.304;
319  }
320  else
321  {
322  aL = 58.83;
323  b1L = 160.7;
324  b2L = -203.5;
325  aT = 0.1452;
326  b1T = 0.02681;
327  b2T = -0.03411;
328  }
329  // laminar flow friction factor for bare Pin bundle - Corner subchannel
330  cL = aL + b1L * (w_over_d - 1) + b2L * std::pow((w_over_d - 1), 2);
331  // turbulent flow friction factor for bare Pin bundle - Corner subchannel
332  cT = aT + b1T * (w_over_d - 1) + b2T * std::pow((w_over_d - 1), 2);
333  }
334 
335  // laminar friction factor
336  auto fL = cL / Re;
337  // turbulent friction factor
338  auto fT = cT / std::pow(Re, 0.18);
339 
340  if (Re < ReL)
341  {
342  // laminar flow
343  return fL;
344  }
345  else if (Re > ReT)
346  {
347  // turbulent flow
348  return fT;
349  }
350  else
351  {
352  // transient flow: psi definition uses a Bulk ReT/ReL number, same for all channels
353  return fL * std::pow((1 - psi), 1.0 / 3.0) * (1 - std::pow(psi, lambda)) +
354  fT * std::pow(psi, 1.0 / 3.0);
355  }
356  }
357 }
358 
359 Real
360 QuadSubChannel1PhaseProblem::computeBeta(unsigned int i_gap, unsigned int iz)
361 {
362  auto beta = _beta;
363  if (!_constant_beta)
364  {
365  const Real & pitch = _subchannel_mesh.getPitch();
366  const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
367  auto chans = _subchannel_mesh.getGapChannels(i_gap);
368  unsigned int i_ch = chans.first;
369  unsigned int j_ch = chans.second;
370  auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
371  auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
372  auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
373  auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
374  auto Si_in = (*_S_flow_soln)(node_in_i);
375  auto Sj_in = (*_S_flow_soln)(node_in_j);
376  auto Si_out = (*_S_flow_soln)(node_out_i);
377  auto Sj_out = (*_S_flow_soln)(node_out_j);
378  // crossflow area between channels i,j (dz*gap_width)
379  auto gap = _subchannel_mesh.getGapWidth(iz, i_gap);
380  auto avg_massflux =
381  0.5 * (((*_mdot_soln)(node_in_i) + (*_mdot_soln)(node_in_j)) / (Si_in + Sj_in) +
382  ((*_mdot_soln)(node_out_i) + (*_mdot_soln)(node_out_j)) / (Si_out + Sj_out));
383  auto S_total = Si_in + Sj_in + Si_out + Sj_out;
384  auto Si = 0.5 * (Si_in + Si_out);
385  auto Sj = 0.5 * (Sj_in + Sj_out);
386  auto w_perim_i = 0.5 * ((*_w_perim_soln)(node_in_i) + (*_w_perim_soln)(node_out_i));
387  auto w_perim_j = 0.5 * ((*_w_perim_soln)(node_in_j) + (*_w_perim_soln)(node_out_j));
388  auto avg_mu = (1 / S_total) * ((*_mu_soln)(node_out_i)*Si_out + (*_mu_soln)(node_in_i)*Si_in +
389  (*_mu_soln)(node_out_j)*Sj_out + (*_mu_soln)(node_in_j)*Sj_in);
390  auto avg_hD = 4.0 * (Si + Sj) / (w_perim_i + w_perim_j);
391  auto Re = avg_massflux * avg_hD / avg_mu;
392  Real gamma = 20.0; // empirical constant
393  Real sf = 1.0; // shape factor
394  Real a = 0.18;
395  Real b = 0.2;
396  auto f = a * std::pow(Re, -b); // Rehme 1992 circular tube friction factor
397  auto k = (1 / S_total) *
398  (_fp->k_from_p_T((*_P_soln)(node_out_i) + _P_out, (*_T_soln)(node_out_i)) * Si_out +
399  _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i)) * Si_in +
400  _fp->k_from_p_T((*_P_soln)(node_out_j) + _P_out, (*_T_soln)(node_out_j)) * Sj_out +
401  _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j)) * Sj_in);
402  auto cp = (1 / S_total) *
403  (_fp->cp_from_p_T((*_P_soln)(node_out_i) + _P_out, (*_T_soln)(node_out_i)) * Si_out +
404  _fp->cp_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i)) * Si_in +
405  _fp->cp_from_p_T((*_P_soln)(node_out_j) + _P_out, (*_T_soln)(node_out_j)) * Sj_out +
406  _fp->cp_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j)) * Sj_in);
407  auto Pr = avg_mu * cp / k; // Prandtl number
408  auto Pr_t = Pr * (Re / gamma) * std::sqrt(f / 8.0); // Turbulent Prandtl number
409  auto delta = pitch; // centroid to centroid distance
410  auto L_x = sf * delta; // axial length scale (gap is the lateral length scale)
411  auto lamda = gap / L_x; // aspect ratio
412  auto a_x = 1.0 - 2.0 * lamda * lamda / libMesh::pi; // velocity coefficient
413  auto z_FP_over_D = (2.0 * L_x / pin_diameter) *
414  (1 + (-0.5 * std::log(lamda) + 0.5 * std::log(4.0) - 0.25) * lamda * lamda);
415  auto Str = 1.0 / (0.822 * (gap / pin_diameter) + 0.144); // Strouhal number (Wu & Trupp 1994)
416  auto freq_factor = 2.0 / std::pow(gamma, 2) * std::sqrt(a / 8.0) * (avg_hD / gap);
417  auto rod_mixing = (1 / Pr_t) * lamda;
418  auto axial_mixing = a_x * z_FP_over_D * Str;
419  // Mixing Stanton number: Stg (eq 25,Kim and Chung (2001), eq 19 (Jeong et. al 2005)
420  beta = freq_factor * (rod_mixing + axial_mixing) * std::pow(Re, -b / 2.0);
421  }
422  mooseAssert(beta > 0, "beta should be positive");
423  return beta;
424 }
425 
426 Real
427 QuadSubChannel1PhaseProblem::computeAddedHeatPin(unsigned int i_ch, unsigned int iz)
428 {
429  // Compute axial location of nodes.
430  auto z2 = _z_grid[iz];
431  auto z1 = _z_grid[iz - 1];
432  auto heated_length = _subchannel_mesh.getHeatedLength();
433  auto unheated_length_entry = _subchannel_mesh.getHeatedLengthEntry();
434  if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
435  MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
436  {
437  // Compute the height of this element.
438  auto dz = z2 - z1;
439  if (_pin_mesh_exist)
440  {
441  auto heat_rate_in = 0.0;
442  auto heat_rate_out = 0.0;
443  for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
444  {
445  auto * node_in = _subchannel_mesh.getPinNode(i_pin, iz - 1);
446  auto * node_out = _subchannel_mesh.getPinNode(i_pin, iz);
447  heat_rate_out += 0.25 * (*_q_prime_soln)(node_out);
448  heat_rate_in += 0.25 * (*_q_prime_soln)(node_in);
449  }
450  return (heat_rate_in + heat_rate_out) * dz / 2.0;
451  }
452  else
453  {
454  auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
455  auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
456  return ((*_q_prime_soln)(node_out) + (*_q_prime_soln)(node_in)) * dz / 2.0;
457  }
458  }
459  else
460  return 0.0;
461 }
462 
463 void
465 {
466  unsigned int last_node = (iblock + 1) * _block_size;
467  unsigned int first_node = iblock * _block_size + 1;
468  if (iblock == 0)
469  {
470  for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
471  {
472  auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
473  auto h_out = _fp->h_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
474  if (h_out < 0)
475  {
476  mooseError(
477  name(), " : Calculation of negative Enthalpy h_out = : ", h_out, " Axial Level= : ", 0);
478  }
479  _h_soln->set(node, h_out);
480  }
481  }
482 
483  if (!_implicit_bool)
484  {
485  for (unsigned int iz = first_node; iz < last_node + 1; iz++)
486  {
487  auto dz = _z_grid[iz] - _z_grid[iz - 1];
488  for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
489  {
490  auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
491  auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
492  auto mdot_in = (*_mdot_soln)(node_in);
493  auto h_in = (*_h_soln)(node_in); // J/kg
494  auto volume = dz * (*_S_flow_soln)(node_in);
495  auto mdot_out = (*_mdot_soln)(node_out);
496  auto h_out = 0.0;
497  Real sumWijh = 0.0;
498  Real sumWijPrimeDhij = 0.0;
499  Real added_enthalpy = computeAddedHeatPin(i_ch, iz);
500  // Calculate sum of crossflow into channel i from channels j around i
501  unsigned int counter = 0;
502  for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
503  {
504  auto chans = _subchannel_mesh.getGapChannels(i_gap);
505  unsigned int ii_ch = chans.first;
506  // i is always the smallest and first index in the mapping
507  unsigned int jj_ch = chans.second;
508  auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
509  auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
510  // Define donor enthalpy
511  Real h_star = 0.0;
512  if (_Wij(i_gap, iz) > 0.0)
513  h_star = (*_h_soln)(node_in_i);
514  else if (_Wij(i_gap, iz) < 0.0)
515  h_star = (*_h_soln)(node_in_j);
516  // take care of the sign by applying the map, use donor cell
517  sumWijh += _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * h_star;
518  sumWijPrimeDhij += _WijPrime(i_gap, iz) * (2 * (*_h_soln)(node_in) -
519  (*_h_soln)(node_in_j) - (*_h_soln)(node_in_i));
520  counter++;
521  }
522  h_out = (mdot_in * h_in - sumWijh - sumWijPrimeDhij + added_enthalpy +
523  _TR * _rho_soln->old(node_out) * _h_soln->old(node_out) * volume / _dt) /
524  (mdot_out + _TR * (*_rho_soln)(node_out)*volume / _dt);
525  if (h_out < 0)
526  {
527  mooseError(name(),
528  " : Calculation of negative Enthalpy h_out = : ",
529  h_out,
530  " Axial Level= : ",
531  iz);
532  }
533  _h_soln->set(node_out, h_out); // J/kg
534  }
535  }
536  }
537  else
538  {
539  LibmeshPetscCall(MatZeroEntries(_hc_time_derivative_mat));
540  LibmeshPetscCall(MatZeroEntries(_hc_advective_derivative_mat));
541  LibmeshPetscCall(MatZeroEntries(_hc_cross_derivative_mat));
542  LibmeshPetscCall(VecZeroEntries(_hc_time_derivative_rhs));
543  LibmeshPetscCall(VecZeroEntries(_hc_advective_derivative_rhs));
544  LibmeshPetscCall(VecZeroEntries(_hc_cross_derivative_rhs));
545  LibmeshPetscCall(VecZeroEntries(_hc_added_heat_rhs));
546  LibmeshPetscCall(MatZeroEntries(_hc_sys_h_mat));
547  LibmeshPetscCall(VecZeroEntries(_hc_sys_h_rhs));
548  for (unsigned int iz = first_node; iz < last_node + 1; iz++)
549  {
550  auto dz = _z_grid[iz] - _z_grid[iz - 1];
551  auto iz_ind = iz - first_node;
552  for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
553  {
554  auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
555  auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
556  auto S_in = (*_S_flow_soln)(node_in);
557  auto S_out = (*_S_flow_soln)(node_out);
558  auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
559  auto volume = dz * S_interp;
560 
561  // interpolation weight coefficient
563 
565  if (iz == first_node)
566  {
567  PetscScalar value_vec_tt =
568  -1.0 * _TR * alpha * (*_rho_soln)(node_in) * (*_h_soln)(node_in)*volume / _dt;
569  PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
570  LibmeshPetscCall(
571  VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
572  }
573  else
574  {
575  PetscInt row_tt = i_ch + _n_channels * iz_ind;
576  PetscInt col_tt = i_ch + _n_channels * (iz_ind - 1);
577  PetscScalar value_tt = _TR * alpha * (*_rho_soln)(node_in)*volume / _dt;
578  LibmeshPetscCall(MatSetValues(
579  _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
580  }
581 
582  // Adding diagonal elements
583  PetscInt row_tt = i_ch + _n_channels * iz_ind;
584  PetscInt col_tt = i_ch + _n_channels * iz_ind;
585  PetscScalar value_tt = _TR * (1.0 - alpha) * (*_rho_soln)(node_out)*volume / _dt;
586  LibmeshPetscCall(MatSetValues(
587  _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
588 
589  // Adding RHS elements
590  PetscScalar rho_old_interp =
591  computeInterpolatedValue(_rho_soln->old(node_out), _rho_soln->old(node_in), 0.5);
592  PetscScalar h_old_interp =
593  computeInterpolatedValue(_h_soln->old(node_out), _h_soln->old(node_in), 0.5);
594  PetscScalar value_vec_tt = _TR * rho_old_interp * h_old_interp * volume / _dt;
595  PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
596  LibmeshPetscCall(
597  VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
598 
600  if (iz == first_node)
601  {
602  PetscScalar value_vec_at = (*_mdot_soln)(node_in) * (*_h_soln)(node_in);
603  PetscInt row_vec_at = i_ch + _n_channels * iz_ind;
604  LibmeshPetscCall(VecSetValues(
605  _hc_advective_derivative_rhs, 1, &row_vec_at, &value_vec_at, ADD_VALUES));
606  }
607  else
608  {
609  PetscInt row_at = i_ch + _n_channels * iz_ind;
610  PetscInt col_at = i_ch + _n_channels * (iz_ind - 1);
611  PetscScalar value_at = -1.0 * (*_mdot_soln)(node_in);
612  LibmeshPetscCall(MatSetValues(
613  _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
614  }
615 
616  // Adding diagonal elements
617  PetscInt row_at = i_ch + _n_channels * iz_ind;
618  PetscInt col_at = i_ch + _n_channels * iz_ind;
619  PetscScalar value_at = (*_mdot_soln)(node_out);
620  LibmeshPetscCall(MatSetValues(
621  _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
622 
624  unsigned int counter = 0;
625  unsigned int cross_index = iz; // iz-1;
626  for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
627  {
628  auto chans = _subchannel_mesh.getGapChannels(i_gap);
629  unsigned int ii_ch = chans.first;
630  unsigned int jj_ch = chans.second;
631  auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
632  auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
633  PetscScalar h_star;
634  // figure out donor axial velocity
635  if (_Wij(i_gap, cross_index) > 0.0)
636  {
637  if (iz == first_node)
638  {
639  h_star = (*_h_soln)(node_in_i);
640  PetscScalar value_vec_ct = -1.0 * alpha *
641  _subchannel_mesh.getCrossflowSign(i_ch, counter) *
642  _Wij(i_gap, cross_index) * h_star;
643  PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
644  LibmeshPetscCall(VecSetValues(
645  _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
646  }
647  else
648  {
649  PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
650  _Wij(i_gap, cross_index);
651  PetscInt row_ct = i_ch + _n_channels * iz_ind;
652  PetscInt col_ct = ii_ch + _n_channels * (iz_ind - 1);
653  LibmeshPetscCall(MatSetValues(
654  _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
655  }
656  PetscScalar value_ct = (1.0 - alpha) *
657  _subchannel_mesh.getCrossflowSign(i_ch, counter) *
658  _Wij(i_gap, cross_index);
659  PetscInt row_ct = i_ch + _n_channels * iz_ind;
660  PetscInt col_ct = ii_ch + _n_channels * iz_ind;
661  LibmeshPetscCall(MatSetValues(
662  _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
663  }
664  else if (_Wij(i_gap, cross_index) < 0.0) // _Wij=0 operations not necessary
665  {
666  if (iz == first_node)
667  {
668  h_star = (*_h_soln)(node_in_j);
669  PetscScalar value_vec_ct = -1.0 * alpha *
670  _subchannel_mesh.getCrossflowSign(i_ch, counter) *
671  _Wij(i_gap, cross_index) * h_star;
672  PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
673  LibmeshPetscCall(VecSetValues(
674  _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
675  }
676  else
677  {
678  PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
679  _Wij(i_gap, cross_index);
680  PetscInt row_ct = i_ch + _n_channels * iz_ind;
681  PetscInt col_ct = jj_ch + _n_channels * (iz_ind - 1);
682  LibmeshPetscCall(MatSetValues(
683  _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
684  }
685  PetscScalar value_ct = (1.0 - alpha) *
686  _subchannel_mesh.getCrossflowSign(i_ch, counter) *
687  _Wij(i_gap, cross_index);
688  PetscInt row_ct = i_ch + _n_channels * iz_ind;
689  PetscInt col_ct = jj_ch + _n_channels * iz_ind;
690  LibmeshPetscCall(MatSetValues(
691  _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
692  }
693 
694  // Turbulent cross flows
695  if (iz == first_node)
696  {
697  PetscScalar value_vec_ct =
698  -2.0 * alpha * (*_h_soln)(node_in)*_WijPrime(i_gap, cross_index);
699  value_vec_ct += alpha * (*_h_soln)(node_in_j)*_WijPrime(i_gap, cross_index);
700  value_vec_ct += alpha * (*_h_soln)(node_in_i)*_WijPrime(i_gap, cross_index);
701  PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
702  LibmeshPetscCall(
703  VecSetValues(_hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
704  }
705  else
706  {
707  PetscScalar value_center_ct = 2.0 * alpha * _WijPrime(i_gap, cross_index);
708  PetscInt row_ct = i_ch + _n_channels * iz_ind;
709  PetscInt col_ct = i_ch + _n_channels * (iz_ind - 1);
710  LibmeshPetscCall(MatSetValues(
711  _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
712 
713  PetscScalar value_left_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
714  row_ct = i_ch + _n_channels * iz_ind;
715  col_ct = jj_ch + _n_channels * (iz_ind - 1);
716  LibmeshPetscCall(MatSetValues(
717  _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
718 
719  PetscScalar value_right_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
720  row_ct = i_ch + _n_channels * iz_ind;
721  col_ct = ii_ch + _n_channels * (iz_ind - 1);
722  LibmeshPetscCall(MatSetValues(
723  _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
724  }
725  PetscScalar value_center_ct = 2.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
726  PetscInt row_ct = i_ch + _n_channels * iz_ind;
727  PetscInt col_ct = i_ch + _n_channels * iz_ind;
728  LibmeshPetscCall(MatSetValues(
729  _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
730 
731  PetscScalar value_left_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
732  row_ct = i_ch + _n_channels * iz_ind;
733  col_ct = jj_ch + _n_channels * iz_ind;
734  LibmeshPetscCall(MatSetValues(
735  _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
736 
737  PetscScalar value_right_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
738  row_ct = i_ch + _n_channels * iz_ind;
739  col_ct = ii_ch + _n_channels * iz_ind;
740  LibmeshPetscCall(MatSetValues(
741  _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
742  counter++;
743  }
744 
746  PetscScalar added_enthalpy = computeAddedHeatPin(i_ch, iz);
747  PetscInt row_vec_ht = i_ch + _n_channels * iz_ind;
748  LibmeshPetscCall(
749  VecSetValues(_hc_added_heat_rhs, 1, &row_vec_ht, &added_enthalpy, ADD_VALUES));
750  }
751  }
753  LibmeshPetscCall(MatAssemblyBegin(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
754  LibmeshPetscCall(MatAssemblyEnd(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
755  LibmeshPetscCall(MatAssemblyBegin(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
756  LibmeshPetscCall(MatAssemblyEnd(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
757  LibmeshPetscCall(MatAssemblyBegin(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
758  LibmeshPetscCall(MatAssemblyEnd(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
759  LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
760  LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
761  // Matrix
762 #if !PETSC_VERSION_LESS_THAN(3, 15, 0)
763  LibmeshPetscCall(MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, UNKNOWN_NONZERO_PATTERN));
764  LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
765  LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
766  LibmeshPetscCall(
767  MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, UNKNOWN_NONZERO_PATTERN));
768  LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
769  LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
770  LibmeshPetscCall(
771  MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, UNKNOWN_NONZERO_PATTERN));
772 #else
773  LibmeshPetscCall(
774  MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
775  LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
776  LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
777  LibmeshPetscCall(
778  MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
779  LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
780  LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
781  LibmeshPetscCall(
782  MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, DIFFERENT_NONZERO_PATTERN));
783 #endif
784  LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
785  LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
787  _console << "Block: " << iblock << " - Enthalpy conservation matrix assembled" << std::endl;
788  // RHS
789  LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_time_derivative_rhs));
790  LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_advective_derivative_rhs));
791  LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_cross_derivative_rhs));
792  LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_added_heat_rhs));
793 
795  {
796  // Assembly the matrix system
797  KSP ksploc;
798  PC pc;
799  Vec sol;
800  LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &sol));
801  LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
802  LibmeshPetscCall(KSPSetOperators(ksploc, _hc_sys_h_mat, _hc_sys_h_mat));
803  LibmeshPetscCall(KSPGetPC(ksploc, &pc));
804  LibmeshPetscCall(PCSetType(pc, PCJACOBI));
805  LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit));
806  LibmeshPetscCall(KSPSetOptionsPrefix(ksploc, "h_sys_"));
807  LibmeshPetscCall(KSPSetFromOptions(ksploc));
808  LibmeshPetscCall(KSPSolve(ksploc, _hc_sys_h_rhs, sol));
809  PetscScalar * xx;
810  LibmeshPetscCall(VecGetArray(sol, &xx));
811  for (unsigned int iz = first_node; iz < last_node + 1; iz++)
812  {
813  auto iz_ind = iz - first_node;
814  for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
815  {
816  auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
817  auto h_out = xx[iz_ind * _n_channels + i_ch];
818  if (h_out < 0)
819  {
820  mooseError(name(),
821  " : Calculation of negative Enthalpy h_out = : ",
822  h_out,
823  " Axial Level= : ",
824  iz);
825  }
826  _h_soln->set(node_out, h_out);
827  }
828  }
829  LibmeshPetscCall(KSPDestroy(&ksploc));
830  LibmeshPetscCall(VecDestroy(&sol));
831  }
832  }
833 }
virtual const Real & getCrossflowSign(unsigned int i_chan, unsigned int i_local) const override
Return a sign for the crossflow given a subchannel index and local neighbor index.
const bool _pin_mesh_exist
Flag that informs if there is a pin mesh or not.
virtual const Real & getPinDiameter() const
Return Pin diameter.
std::unique_ptr< SolutionHandle > _T_soln
virtual Node * getPinNode(unsigned int i_pin, unsigned iz) const override
Get the pin mesh node for a given pin index and elevation index.
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.
T & getMesh(MooseMesh &mesh)
function to cast mesh
Definition: SCM.h:35
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
Real _TR
Flag that activates or deactivates the transient parts of the equations we solve by multiplication...
const PostprocessorValue & _P_out
Outlet Pressure.
virtual Real computeAddedHeatPin(unsigned int i_ch, unsigned int iz) override
Computes added heat for channel i_ch and cell iz.
const bool _constant_beta
Flag that activates the use of constant beta.
static InputParameters validParams()
virtual const std::vector< unsigned int > & getChannelGaps(unsigned int i_chan) const override
Return a vector of gap indices for a given channel index.
std::unique_ptr< SolutionHandle > _S_flow_soln
virtual const std::pair< unsigned int, unsigned int > & getGapPins(unsigned int i_gap) const override
Return a pair of pin indices for a given gap index.
int delta(unsigned int i, unsigned int j)
Delta function, which returns zero if $i j$ and unity if $i=j$.
virtual void computeh(int iblock) override
Computes Enthalpy per channel for block iblock.
static InputParameters validParams()
Creates the mesh of subchannels in a quadrilateral lattice.
PetscScalar computeInterpolatedValue(PetscScalar topValue, PetscScalar botValue, PetscScalar Peclet=0.0)
const bool _segregated_bool
Segregated solve.
virtual const std::string & name() const
void mooseWarning(Args &&... args) const
void addRequiredParam(const std::string &name, const std::string &doc_string)
virtual Real getGapWidth(unsigned int axial_index, unsigned int gap_index) const override
Return gap width for a given gap index.
virtual const std::pair< unsigned int, unsigned int > & getGapChannels(unsigned int i_gap) const override
Return a pair of subchannel indices for a given gap index.
registerMooseObject("SubChannelApp", QuadSubChannel1PhaseProblem)
std::unique_ptr< SolutionHandle > _rho_soln
const Real & getGap() const
Returns the gap, not to be confused with the gap between pins, this refers to the gap next to the duc...
const bool _default_friction_model
Flag that activates one of the two friction models (default: f=a*Re^b, non-default: Todreas-Kazimi) ...
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
virtual const Real & getHeatedLength() const
Return heated length.
Real f(Real x)
Test function for Brents method.
std::unique_ptr< SolutionHandle > _q_prime_soln
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.
virtual Real computeFrictionFactor(FrictionStruct friction_args) override
Returns friction factor.
std::unique_ptr< SolutionHandle > _Dpin_soln
virtual Real computeBeta(unsigned int i_gap, unsigned int iz) override
Computes turbulent mixing coefficient.
virtual const std::vector< Real > & getReductionBlockage() const
Get area reduction of blocked subchannels.
Real volume(const MeshBase &mesh, unsigned int dim=libMesh::invalid_uint)
Quadrilateral subchannel solver.
const Real & _beta
Thermal diffusion coefficient used in turbulent crossflow.
Base class for the 1-phase steady-state/transient subchannel solver.
std::unique_ptr< SolutionHandle > _mdot_soln
virtual const Real & getPitch() const override
Return the pitch between 2 subchannels.
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
std::unique_ptr< SolutionHandle > _displacement_soln
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 Node * getChannelNode(unsigned int i_chan, unsigned iz) const override
Get the subchannel mesh node for a given channel index and elevation index.
std::vector< std::vector< Real > > _gij_map
Vector to store gap size.
virtual const std::vector< unsigned int > & getPinChannels(unsigned int i_pin) const override
Return a vector of channel indices for a given Pin index.
const bool _verbose_subchannel
Boolean to printout information related to subchannel solve.
void mooseError(Args &&... args) const
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
virtual EChannelType getSubchannelType(unsigned int index) const override
Return the type of the subchannel for given subchannel index.
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
virtual const std::vector< unsigned int > & getChannelPins(unsigned int i_chan) const override
Return a vector of pin indices for a given channel index.
MooseUnits pow(const MooseUnits &, int)
QuadSubChannel1PhaseProblem(const InputParameters &params)
std::unique_ptr< SolutionHandle > _P_soln
virtual void initializeSolution() override
Function to initialize the solution & geometry fields.
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
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