Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "QuadSubChannel1PhaseProblem.h"
11 : #include "AuxiliarySystem.h"
12 : #include "SCM.h"
13 :
14 : registerMooseObject("SubChannelApp", QuadSubChannel1PhaseProblem);
15 :
16 : InputParameters
17 446 : QuadSubChannel1PhaseProblem::validParams()
18 : {
19 446 : InputParameters params = SubChannel1PhaseProblem::validParams();
20 446 : params.addClassDescription(
21 : "Solver class for subchannels in a square lattice assembly and bare fuel pins");
22 892 : params.addRequiredParam<Real>("beta",
23 : "Thermal diffusion coefficient used in turbulent crossflow.");
24 892 : params.addParam<bool>(
25 : "default_friction_model",
26 892 : 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 892 : params.addParam<bool>(
30 : "constant_beta",
31 892 : true,
32 : "Boolean to define the use of a constant beta or beta correlation (Kim and Chung, 2001)");
33 446 : return params;
34 0 : }
35 :
36 223 : QuadSubChannel1PhaseProblem::QuadSubChannel1PhaseProblem(const InputParameters & params)
37 : : SubChannel1PhaseProblem(params),
38 446 : _subchannel_mesh(SCM::getMesh<QuadSubChannelMesh>(_mesh)),
39 446 : _beta(getParam<Real>("beta")),
40 446 : _default_friction_model(getParam<bool>("default_friction_model")),
41 669 : _constant_beta(getParam<bool>("constant_beta"))
42 : {
43 223 : }
44 :
45 : void
46 297 : QuadSubChannel1PhaseProblem::initializeSolution()
47 : {
48 : /// update surface area, wetted perimeter based on: Dpin, displacement
49 297 : if (_deformation)
50 : {
51 : Real standard_area, additional_area, wetted_perimeter, displaced_area;
52 9 : auto pitch = _subchannel_mesh.getPitch();
53 9 : auto side_gap = _subchannel_mesh.getSideGap();
54 9 : auto z_blockage = _subchannel_mesh.getZBlockage();
55 9 : auto index_blockage = _subchannel_mesh.getIndexBlockage();
56 9 : auto reduction_blockage = _subchannel_mesh.getReductionBlockage();
57 138 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
58 : {
59 4665 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
60 : {
61 4536 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
62 4536 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
63 4536 : auto Z = _z_grid[iz];
64 : Real rod_area = 0.0;
65 : Real rod_perimeter = 0.0;
66 17136 : for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
67 : {
68 12600 : auto * pin_node = _subchannel_mesh.getPinNode(i_pin, iz);
69 12600 : rod_area += 0.25 * 0.25 * M_PI * (*_Dpin_soln)(pin_node) * (*_Dpin_soln)(pin_node);
70 12600 : rod_perimeter += 0.25 * M_PI * (*_Dpin_soln)(pin_node);
71 : }
72 :
73 4536 : if (subch_type == EChannelType::CORNER)
74 : {
75 504 : standard_area = 0.25 * pitch * pitch;
76 504 : displaced_area = (2 * side_gap + pitch) * (*_displacement_soln)(node) / sqrt(2) +
77 504 : (*_displacement_soln)(node) * (*_displacement_soln)(node) / 2;
78 504 : additional_area = pitch * side_gap + side_gap * side_gap;
79 504 : wetted_perimeter =
80 504 : rod_perimeter + pitch + 2 * side_gap + 2 * (*_displacement_soln)(node) / sqrt(2);
81 : }
82 4032 : else if (subch_type == EChannelType::EDGE)
83 : {
84 2016 : standard_area = 0.5 * pitch * pitch;
85 2016 : additional_area = pitch * side_gap;
86 2016 : displaced_area = pitch * (*_displacement_soln)(node);
87 2016 : wetted_perimeter = rod_perimeter + pitch;
88 : }
89 : else
90 : {
91 2016 : standard_area = pitch * pitch;
92 : displaced_area = 0.0;
93 : additional_area = 0.0;
94 : wetted_perimeter = rod_perimeter;
95 : }
96 :
97 : /// Calculate subchannel area
98 4536 : auto subchannel_area = displaced_area + standard_area + additional_area - rod_area;
99 :
100 : /// Correct subchannel area and wetted perimeter in case of overlapping pins
101 : auto overlapping_pin_area = 0.0;
102 : auto overlapping_wetted_perimeter = 0.0;
103 19656 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
104 : {
105 15120 : auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
106 : auto pin_1 = gap_pins.first;
107 : auto pin_2 = gap_pins.second;
108 15120 : auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
109 15120 : auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
110 15120 : auto Diameter1 = (*_Dpin_soln)(pin_node_1);
111 15120 : auto Radius1 = Diameter1 / 2.0;
112 15120 : auto Diameter2 = (*_Dpin_soln)(pin_node_2);
113 15120 : auto Radius2 = Diameter2 / 2.0;
114 15120 : auto pitch = _subchannel_mesh.getPitch();
115 :
116 15120 : if (pitch < (Radius1 + Radius2)) // overlapping pins
117 : {
118 0 : mooseWarning(" The gap of index : '", i_gap, " at axial cell ", iz, " ' is blocked.");
119 0 : auto cos1 =
120 0 : (pitch * pitch + Radius1 * Radius1 - Radius2 * Radius2) / (2 * pitch * Radius1);
121 0 : auto cos2 =
122 0 : (pitch * pitch + Radius2 * Radius2 - Radius1 * Radius1) / (2 * pitch * Radius2);
123 0 : auto angle1 = 2.0 * acos(cos1);
124 0 : auto angle2 = 2.0 * acos(cos2);
125 : // half of the intersecting arc-length
126 0 : overlapping_wetted_perimeter += 0.5 * angle1 * Radius1 + 0.5 * angle2 * Radius2;
127 : // Half of the overlapping area
128 0 : overlapping_pin_area +=
129 0 : 0.5 * Radius1 * Radius1 * acos(cos1) + 0.5 * Radius2 * Radius2 * acos(cos2) -
130 0 : 0.25 * sqrt((-pitch + Radius1 + Radius2) * (pitch + Radius1 - Radius2) *
131 0 : (pitch - Radius1 + Radius2) * (pitch + Radius1 + Radius2));
132 : }
133 : }
134 4536 : subchannel_area += overlapping_pin_area; // correct surface area
135 4536 : wetted_perimeter += -overlapping_wetted_perimeter; // correct wetted perimeter
136 :
137 : /// Apply area reduction on subchannels affected by blockage
138 : auto index = 0;
139 9072 : for (const auto & i_blockage : index_blockage)
140 : {
141 4536 : if (i_ch == i_blockage && (Z >= z_blockage.front() && Z <= z_blockage.back()))
142 : {
143 6 : subchannel_area *= reduction_blockage[index];
144 : }
145 4536 : index++;
146 : }
147 :
148 4536 : _S_flow_soln->set(node, subchannel_area);
149 4536 : _w_perim_soln->set(node, wetted_perimeter);
150 : }
151 : }
152 : /// update map of gap between pins (gij) based on: Dpin, displacement
153 138 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
154 : {
155 7689 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
156 : {
157 7560 : auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
158 : auto pin_1 = gap_pins.first;
159 : auto pin_2 = gap_pins.second;
160 7560 : auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
161 7560 : auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
162 7560 : if (pin_1 == pin_2) // Corner or edge gap
163 : {
164 : auto displacement = 0.0;
165 : auto counter = 0.0;
166 12600 : for (auto i_ch : _subchannel_mesh.getPinChannels(pin_1))
167 : {
168 10080 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
169 10080 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
170 10080 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
171 : {
172 6048 : displacement += (*_displacement_soln)(node);
173 6048 : counter += 1.0;
174 : }
175 : }
176 2520 : displacement = displacement / counter;
177 2520 : _subchannel_mesh._gij_map[iz][i_gap] =
178 2520 : (pitch - (*_Dpin_soln)(pin_node_1)) / 2.0 + side_gap + displacement;
179 : }
180 : else // center gap
181 : {
182 5040 : _subchannel_mesh._gij_map[iz][i_gap] =
183 5040 : 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 7560 : if (_subchannel_mesh._gij_map[iz][i_gap] <= 0.0)
187 0 : _subchannel_mesh._gij_map[iz][i_gap] = 0.0;
188 : }
189 : }
190 9 : }
191 :
192 7647 : for (unsigned int iz = 1; iz < _n_cells + 1; iz++)
193 : {
194 120656 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
195 : {
196 113306 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
197 113306 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
198 113306 : _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 297 : _aux->solution().close();
205 297 : }
206 :
207 : Real
208 13205844 : QuadSubChannel1PhaseProblem::computeFrictionFactor(FrictionStruct friction_args)
209 : {
210 13205844 : auto Re = friction_args.Re;
211 13205844 : auto i_ch = friction_args.i_ch;
212 : /// Pang, B. et al. KIT, 2013
213 13205844 : if (_default_friction_model)
214 : {
215 : Real a, b;
216 10784724 : if (Re < 1)
217 : {
218 : return 64.0;
219 : }
220 10784724 : else if (Re >= 1 and Re < 5000)
221 : {
222 : a = 64.0;
223 : b = -1.0;
224 : }
225 8077284 : 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 10784724 : return a * std::pow(Re, b);
236 : }
237 : /// Todreas-Kazimi NUCLEAR SYSTEMS, second edition, Volume 1, 2011
238 : else
239 : {
240 : Real aL, b1L, b2L, cL;
241 : Real aT, b1T, b2T, cT;
242 2421120 : auto pitch = _subchannel_mesh.getPitch();
243 2421120 : 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 2421120 : auto side_gap = _subchannel_mesh.getSideGap();
247 2421120 : auto w = (pin_diameter / 2.0) + (pitch / 2.0) + side_gap;
248 2421120 : auto p_over_d = pitch / pin_diameter;
249 2421120 : auto w_over_d = w / pin_diameter;
250 2421120 : auto ReL = std::pow(10, (p_over_d - 1)) * 320.0;
251 2421120 : auto ReT = std::pow(10, 0.7 * (p_over_d - 1)) * 1.0E+4;
252 2421120 : auto psi = std::log(Re / ReL) / std::log(ReT / ReL);
253 2421120 : 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 2421120 : if (subch_type == EChannelType::CENTER)
260 : {
261 1116720 : 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 1116720 : 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 1116720 : cT = aT + b1T * (p_over_d - 1) + b2T * std::pow((p_over_d - 1), 2);
283 : }
284 1304400 : else if (subch_type == EChannelType::EDGE)
285 : {
286 1043520 : 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 1043520 : 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 1043520 : cT = aT + b1T * (w_over_d - 1) + b2T * std::pow((w_over_d - 1), 2);
308 : }
309 : else
310 : {
311 260880 : 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 260880 : 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 260880 : cT = aT + b1T * (w_over_d - 1) + b2T * std::pow((w_over_d - 1), 2);
333 : }
334 :
335 : // laminar friction factor
336 2421120 : auto fL = cL / Re;
337 : // turbulent friction factor
338 2421120 : auto fT = cT / std::pow(Re, 0.18);
339 :
340 2421120 : if (Re < ReL)
341 : {
342 : // laminar flow
343 : return fL;
344 : }
345 2421120 : 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 0 : return fL * std::pow((1 - psi), 1.0 / 3.0) * (1 - std::pow(psi, lambda)) +
354 0 : fT * std::pow(psi, 1.0 / 3.0);
355 : }
356 : }
357 : }
358 :
359 : Real
360 21575664 : QuadSubChannel1PhaseProblem::computeBeta(unsigned int i_gap, unsigned int iz, bool /*enthalpy*/)
361 : {
362 21575664 : auto beta = _beta;
363 21575664 : if (!_constant_beta)
364 : {
365 3913200 : const Real & pitch = _subchannel_mesh.getPitch();
366 3913200 : const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
367 3913200 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
368 : unsigned int i_ch = chans.first;
369 : unsigned int j_ch = chans.second;
370 3913200 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
371 3913200 : auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
372 3913200 : auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
373 3913200 : auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
374 3913200 : auto Si_in = (*_S_flow_soln)(node_in_i);
375 3913200 : auto Sj_in = (*_S_flow_soln)(node_in_j);
376 3913200 : auto Si_out = (*_S_flow_soln)(node_out_i);
377 3913200 : auto Sj_out = (*_S_flow_soln)(node_out_j);
378 : // crossflow area between channels i,j (dz*gap_width)
379 3913200 : auto gap = _subchannel_mesh.getGapWidth(iz, i_gap);
380 : auto avg_massflux =
381 3913200 : 0.5 * (((*_mdot_soln)(node_in_i) + (*_mdot_soln)(node_in_j)) / (Si_in + Sj_in) +
382 3913200 : ((*_mdot_soln)(node_out_i) + (*_mdot_soln)(node_out_j)) / (Si_out + Sj_out));
383 3913200 : auto S_total = Si_in + Sj_in + Si_out + Sj_out;
384 3913200 : auto Si = 0.5 * (Si_in + Si_out);
385 3913200 : auto Sj = 0.5 * (Sj_in + Sj_out);
386 3913200 : auto w_perim_i = 0.5 * ((*_w_perim_soln)(node_in_i) + (*_w_perim_soln)(node_out_i));
387 3913200 : auto w_perim_j = 0.5 * ((*_w_perim_soln)(node_in_j) + (*_w_perim_soln)(node_out_j));
388 3913200 : auto avg_mu = (1 / S_total) * ((*_mu_soln)(node_out_i)*Si_out + (*_mu_soln)(node_in_i)*Si_in +
389 3913200 : (*_mu_soln)(node_out_j)*Sj_out + (*_mu_soln)(node_in_j)*Sj_in);
390 3913200 : auto avg_hD = 4.0 * (Si + Sj) / (w_perim_i + w_perim_j);
391 3913200 : 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 3913200 : auto f = a * std::pow(Re, -b); // Rehme 1992 circular tube friction factor
397 : auto k = (1 / S_total) *
398 3913200 : (_fp->k_from_p_T((*_P_soln)(node_out_i) + _P_out, (*_T_soln)(node_out_i)) * Si_out +
399 3913200 : _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i)) * Si_in +
400 3913200 : _fp->k_from_p_T((*_P_soln)(node_out_j) + _P_out, (*_T_soln)(node_out_j)) * Sj_out +
401 3913200 : _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 3913200 : (_fp->cp_from_p_T((*_P_soln)(node_out_i) + _P_out, (*_T_soln)(node_out_i)) * Si_out +
404 3913200 : _fp->cp_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i)) * Si_in +
405 3913200 : _fp->cp_from_p_T((*_P_soln)(node_out_j) + _P_out, (*_T_soln)(node_out_j)) * Sj_out +
406 3913200 : _fp->cp_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j)) * Sj_in);
407 3913200 : auto Pr = avg_mu * cp / k; // Prandtl number
408 3913200 : auto Pr_t = Pr * (Re / gamma) * std::sqrt(f / 8.0); // Turbulent Prandtl number
409 3913200 : auto delta = pitch; // centroid to centroid distance
410 : auto L_x = sf * delta; // axial length scale (gap is the lateral length scale)
411 3913200 : auto lamda = gap / L_x; // aspect ratio
412 3913200 : auto a_x = 1.0 - 2.0 * lamda * lamda / libMesh::pi; // velocity coefficient
413 3913200 : auto z_FP_over_D = (2.0 * L_x / pin_diameter) *
414 3913200 : (1 + (-0.5 * std::log(lamda) + 0.5 * std::log(4.0) - 0.25) * lamda * lamda);
415 3913200 : auto Str = 1.0 / (0.822 * (gap / pin_diameter) + 0.144); // Strouhal number (Wu & Trupp 1994)
416 3913200 : auto freq_factor = 2.0 / std::pow(gamma, 2) * std::sqrt(a / 8.0) * (avg_hD / gap);
417 3913200 : auto rod_mixing = (1 / Pr_t) * lamda;
418 3913200 : 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 3913200 : beta = freq_factor * (rod_mixing + axial_mixing) * std::pow(Re, -b / 2.0);
421 : }
422 : mooseAssert(beta > 0, "beta should be positive");
423 21575664 : return beta;
424 : }
425 :
426 : Real
427 332592 : QuadSubChannel1PhaseProblem::computeAddedHeatPin(unsigned int i_ch, unsigned int iz)
428 : {
429 : // Compute axial location of nodes.
430 332592 : auto z2 = _z_grid[iz];
431 332592 : auto z1 = _z_grid[iz - 1];
432 332592 : auto heated_length = _subchannel_mesh.getHeatedLength();
433 332592 : auto unheated_length_entry = _subchannel_mesh.getHeatedLengthEntry();
434 332592 : if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
435 328326 : MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
436 : {
437 : // Compute the height of this element.
438 324060 : auto dz = z2 - z1;
439 324060 : if (_pin_mesh_exist)
440 : {
441 : auto heat_rate_in = 0.0;
442 : auto heat_rate_out = 0.0;
443 1053360 : for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
444 : {
445 747840 : auto * node_in = _subchannel_mesh.getPinNode(i_pin, iz - 1);
446 747840 : auto * node_out = _subchannel_mesh.getPinNode(i_pin, iz);
447 747840 : heat_rate_out += 0.25 * (*_q_prime_soln)(node_out);
448 747840 : heat_rate_in += 0.25 * (*_q_prime_soln)(node_in);
449 : }
450 305520 : return (heat_rate_in + heat_rate_out) * dz / 2.0;
451 : }
452 : else
453 : {
454 18540 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
455 18540 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
456 18540 : 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 : Real
464 0 : QuadSubChannel1PhaseProblem::getSubChannelPeripheralDuctWidth(unsigned int i_ch)
465 : {
466 0 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
467 0 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
468 : {
469 0 : auto width = _subchannel_mesh.getPitch();
470 0 : if (subch_type == EChannelType::CORNER)
471 0 : width += 2.0 * _subchannel_mesh.getSideGap();
472 0 : return width;
473 : }
474 : else
475 0 : mooseError("Channel is not a perimetric subchannel ");
476 : }
477 :
478 : void
479 722 : QuadSubChannel1PhaseProblem::computeh(int iblock)
480 : {
481 722 : unsigned int last_node = (iblock + 1) * _block_size;
482 722 : unsigned int first_node = iblock * _block_size + 1;
483 722 : if (iblock == 0)
484 : {
485 16064 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
486 : {
487 15342 : auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
488 15342 : auto h_out = _fp->h_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
489 15342 : if (h_out < 0)
490 : {
491 0 : mooseError(
492 0 : name(), " : Calculation of negative Enthalpy h_out = : ", h_out, " Axial Level= : ", 0);
493 : }
494 15342 : _h_soln->set(node, h_out);
495 : }
496 : }
497 :
498 722 : if (!_implicit_bool)
499 : {
500 8308 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
501 : {
502 8060 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
503 50780 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
504 : {
505 42720 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
506 42720 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
507 42720 : auto mdot_in = (*_mdot_soln)(node_in);
508 42720 : auto h_in = (*_h_soln)(node_in); // J/kg
509 42720 : auto volume = dz * (*_S_flow_soln)(node_in);
510 42720 : auto mdot_out = (*_mdot_soln)(node_out);
511 42720 : auto h_out = 0.0;
512 : Real sumWijh = 0.0;
513 : Real sumWijPrimeDhij = 0.0;
514 42720 : Real added_enthalpy = computeAddedHeatPin(i_ch, iz);
515 : // Calculate sum of crossflow into channel i from channels j around i
516 : unsigned int counter = 0;
517 148080 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
518 : {
519 105360 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
520 : unsigned int ii_ch = chans.first;
521 : // i is always the smallest and first index in the mapping
522 : unsigned int jj_ch = chans.second;
523 105360 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
524 105360 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
525 : // Define donor enthalpy
526 : Real h_star = 0.0;
527 105360 : if (_Wij(i_gap, iz) > 0.0)
528 42408 : h_star = (*_h_soln)(node_in_i);
529 62952 : else if (_Wij(i_gap, iz) < 0.0)
530 35352 : h_star = (*_h_soln)(node_in_j);
531 : // take care of the sign by applying the map, use donor cell
532 105360 : sumWijh += _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * h_star;
533 105360 : sumWijPrimeDhij += _WijPrime(i_gap, iz) * (2 * (*_h_soln)(node_in) -
534 105360 : (*_h_soln)(node_in_j) - (*_h_soln)(node_in_i));
535 105360 : counter++;
536 : }
537 128160 : h_out = (mdot_in * h_in - sumWijh - sumWijPrimeDhij + added_enthalpy +
538 42720 : _TR * _rho_soln->old(node_out) * _h_soln->old(node_out) * volume / _dt) /
539 42720 : (mdot_out + _TR * (*_rho_soln)(node_out)*volume / _dt);
540 42720 : if (h_out < 0)
541 : {
542 0 : mooseError(name(),
543 : " : Calculation of negative Enthalpy h_out = : ",
544 : h_out,
545 : " Axial Level= : ",
546 : iz);
547 : }
548 42720 : _h_soln->set(node_out, h_out); // J/kg
549 : }
550 : }
551 : }
552 : else
553 : {
554 474 : LibmeshPetscCall(MatZeroEntries(_hc_time_derivative_mat));
555 474 : LibmeshPetscCall(MatZeroEntries(_hc_advective_derivative_mat));
556 474 : LibmeshPetscCall(MatZeroEntries(_hc_cross_derivative_mat));
557 474 : LibmeshPetscCall(VecZeroEntries(_hc_time_derivative_rhs));
558 474 : LibmeshPetscCall(VecZeroEntries(_hc_advective_derivative_rhs));
559 474 : LibmeshPetscCall(VecZeroEntries(_hc_cross_derivative_rhs));
560 474 : LibmeshPetscCall(VecZeroEntries(_hc_added_heat_rhs));
561 474 : LibmeshPetscCall(MatZeroEntries(_hc_sys_h_mat));
562 474 : LibmeshPetscCall(VecZeroEntries(_hc_sys_h_rhs));
563 12750 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
564 : {
565 12276 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
566 12276 : auto iz_ind = iz - first_node;
567 299340 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
568 : {
569 287064 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
570 287064 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
571 287064 : auto S_in = (*_S_flow_soln)(node_in);
572 287064 : auto S_out = (*_S_flow_soln)(node_out);
573 287064 : auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
574 287064 : auto volume = dz * S_interp;
575 :
576 : // interpolation weight coefficient
577 287064 : auto alpha = computeInterpolationCoefficients(0.5);
578 :
579 : /// Time derivative term
580 287064 : if (iz == first_node)
581 : {
582 : PetscScalar value_vec_tt =
583 12420 : -1.0 * _TR * alpha * (*_rho_soln)(node_in) * (*_h_soln)(node_in)*volume / _dt;
584 12420 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
585 12420 : LibmeshPetscCall(
586 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
587 : }
588 : else
589 : {
590 274644 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
591 274644 : PetscInt col_tt = i_ch + _n_channels * (iz_ind - 1);
592 274644 : PetscScalar value_tt = _TR * alpha * (*_rho_soln)(node_in)*volume / _dt;
593 274644 : LibmeshPetscCall(MatSetValues(
594 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
595 : }
596 :
597 : // Adding diagonal elements
598 287064 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
599 287064 : PetscInt col_tt = i_ch + _n_channels * iz_ind;
600 287064 : PetscScalar value_tt = _TR * (1.0 - alpha) * (*_rho_soln)(node_out)*volume / _dt;
601 287064 : LibmeshPetscCall(MatSetValues(
602 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
603 :
604 : // Adding RHS elements
605 : PetscScalar rho_old_interp =
606 287064 : computeInterpolatedValue(_rho_soln->old(node_out), _rho_soln->old(node_in), 0.5);
607 : PetscScalar h_old_interp =
608 287064 : computeInterpolatedValue(_h_soln->old(node_out), _h_soln->old(node_in), 0.5);
609 287064 : PetscScalar value_vec_tt = _TR * rho_old_interp * h_old_interp * volume / _dt;
610 287064 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
611 287064 : LibmeshPetscCall(
612 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
613 :
614 : /// Advective derivative term
615 287064 : if (iz == first_node)
616 : {
617 12420 : PetscScalar value_vec_at = (*_mdot_soln)(node_in) * (*_h_soln)(node_in);
618 12420 : PetscInt row_vec_at = i_ch + _n_channels * iz_ind;
619 12420 : LibmeshPetscCall(VecSetValues(
620 : _hc_advective_derivative_rhs, 1, &row_vec_at, &value_vec_at, ADD_VALUES));
621 : }
622 : else
623 : {
624 274644 : PetscInt row_at = i_ch + _n_channels * iz_ind;
625 274644 : PetscInt col_at = i_ch + _n_channels * (iz_ind - 1);
626 274644 : PetscScalar value_at = -1.0 * (*_mdot_soln)(node_in);
627 274644 : LibmeshPetscCall(MatSetValues(
628 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
629 : }
630 :
631 : // Adding diagonal elements
632 287064 : PetscInt row_at = i_ch + _n_channels * iz_ind;
633 287064 : PetscInt col_at = i_ch + _n_channels * iz_ind;
634 287064 : PetscScalar value_at = (*_mdot_soln)(node_out);
635 287064 : LibmeshPetscCall(MatSetValues(
636 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
637 :
638 : /// Cross derivative term
639 : unsigned int counter = 0;
640 : unsigned int cross_index = iz; // iz-1;
641 1186488 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
642 : {
643 899424 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
644 : unsigned int ii_ch = chans.first;
645 : unsigned int jj_ch = chans.second;
646 899424 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
647 899424 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
648 : PetscScalar h_star;
649 : // figure out donor axial velocity
650 899424 : if (_Wij(i_gap, cross_index) > 0.0)
651 : {
652 529524 : if (iz == first_node)
653 : {
654 22320 : h_star = (*_h_soln)(node_in_i);
655 44640 : PetscScalar value_vec_ct = -1.0 * alpha *
656 22320 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
657 22320 : _Wij(i_gap, cross_index) * h_star;
658 22320 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
659 22320 : LibmeshPetscCall(VecSetValues(
660 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
661 : }
662 : else
663 : {
664 507204 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
665 507204 : _Wij(i_gap, cross_index);
666 507204 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
667 507204 : PetscInt col_ct = ii_ch + _n_channels * (iz_ind - 1);
668 507204 : LibmeshPetscCall(MatSetValues(
669 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
670 : }
671 529524 : PetscScalar value_ct = (1.0 - alpha) *
672 529524 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
673 529524 : _Wij(i_gap, cross_index);
674 529524 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
675 529524 : PetscInt col_ct = ii_ch + _n_channels * iz_ind;
676 529524 : LibmeshPetscCall(MatSetValues(
677 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
678 : }
679 369900 : else if (_Wij(i_gap, cross_index) < 0.0) // _Wij=0 operations not necessary
680 : {
681 319212 : if (iz == first_node)
682 : {
683 16128 : h_star = (*_h_soln)(node_in_j);
684 32256 : PetscScalar value_vec_ct = -1.0 * alpha *
685 16128 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
686 16128 : _Wij(i_gap, cross_index) * h_star;
687 16128 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
688 16128 : LibmeshPetscCall(VecSetValues(
689 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
690 : }
691 : else
692 : {
693 303084 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
694 303084 : _Wij(i_gap, cross_index);
695 303084 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
696 303084 : PetscInt col_ct = jj_ch + _n_channels * (iz_ind - 1);
697 303084 : LibmeshPetscCall(MatSetValues(
698 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
699 : }
700 319212 : PetscScalar value_ct = (1.0 - alpha) *
701 319212 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
702 319212 : _Wij(i_gap, cross_index);
703 319212 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
704 319212 : PetscInt col_ct = jj_ch + _n_channels * iz_ind;
705 319212 : LibmeshPetscCall(MatSetValues(
706 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
707 : }
708 :
709 : // Turbulent cross flows
710 899424 : if (iz == first_node)
711 : {
712 : PetscScalar value_vec_ct =
713 39888 : -2.0 * alpha * (*_h_soln)(node_in)*_WijPrime(i_gap, cross_index);
714 39888 : value_vec_ct += alpha * (*_h_soln)(node_in_j)*_WijPrime(i_gap, cross_index);
715 39888 : value_vec_ct += alpha * (*_h_soln)(node_in_i)*_WijPrime(i_gap, cross_index);
716 39888 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
717 39888 : LibmeshPetscCall(
718 : VecSetValues(_hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
719 : }
720 : else
721 : {
722 859536 : PetscScalar value_center_ct = 2.0 * alpha * _WijPrime(i_gap, cross_index);
723 859536 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
724 859536 : PetscInt col_ct = i_ch + _n_channels * (iz_ind - 1);
725 859536 : LibmeshPetscCall(MatSetValues(
726 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
727 :
728 859536 : PetscScalar value_left_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
729 859536 : row_ct = i_ch + _n_channels * iz_ind;
730 859536 : col_ct = jj_ch + _n_channels * (iz_ind - 1);
731 859536 : LibmeshPetscCall(MatSetValues(
732 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
733 :
734 859536 : PetscScalar value_right_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
735 859536 : row_ct = i_ch + _n_channels * iz_ind;
736 859536 : col_ct = ii_ch + _n_channels * (iz_ind - 1);
737 859536 : LibmeshPetscCall(MatSetValues(
738 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
739 : }
740 899424 : PetscScalar value_center_ct = 2.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
741 899424 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
742 899424 : PetscInt col_ct = i_ch + _n_channels * iz_ind;
743 899424 : LibmeshPetscCall(MatSetValues(
744 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
745 :
746 899424 : PetscScalar value_left_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
747 899424 : row_ct = i_ch + _n_channels * iz_ind;
748 899424 : col_ct = jj_ch + _n_channels * iz_ind;
749 899424 : LibmeshPetscCall(MatSetValues(
750 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
751 :
752 899424 : PetscScalar value_right_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
753 899424 : row_ct = i_ch + _n_channels * iz_ind;
754 899424 : col_ct = ii_ch + _n_channels * iz_ind;
755 899424 : LibmeshPetscCall(MatSetValues(
756 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
757 899424 : counter++;
758 : }
759 :
760 : /// Added heat enthalpy
761 287064 : PetscScalar added_enthalpy = computeAddedHeatPin(i_ch, iz);
762 287064 : PetscInt row_vec_ht = i_ch + _n_channels * iz_ind;
763 287064 : LibmeshPetscCall(
764 : VecSetValues(_hc_added_heat_rhs, 1, &row_vec_ht, &added_enthalpy, ADD_VALUES));
765 : }
766 : }
767 : /// Assembling system
768 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
769 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
770 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
771 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
772 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
773 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
774 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
775 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
776 : // Matrix
777 : #if !PETSC_VERSION_LESS_THAN(3, 15, 0)
778 474 : LibmeshPetscCall(MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, UNKNOWN_NONZERO_PATTERN));
779 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
780 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
781 474 : LibmeshPetscCall(
782 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, UNKNOWN_NONZERO_PATTERN));
783 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
784 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
785 474 : LibmeshPetscCall(
786 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, UNKNOWN_NONZERO_PATTERN));
787 : #else
788 : LibmeshPetscCall(
789 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
790 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
791 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
792 : LibmeshPetscCall(
793 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
794 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
795 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
796 : LibmeshPetscCall(
797 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, DIFFERENT_NONZERO_PATTERN));
798 : #endif
799 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
800 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
801 474 : if (_verbose_subchannel)
802 294 : _console << "Block: " << iblock << " - Enthalpy conservation matrix assembled" << std::endl;
803 : // RHS
804 474 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_time_derivative_rhs));
805 474 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_advective_derivative_rhs));
806 474 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_cross_derivative_rhs));
807 474 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_added_heat_rhs));
808 :
809 474 : if (_segregated_bool || (!_monolithic_thermal_bool))
810 : {
811 : // Assembly the matrix system
812 : KSP ksploc;
813 : PC pc;
814 : Vec sol;
815 360 : LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &sol));
816 360 : LibmeshPetscCall(KSPCreate(PETSC_COMM_SELF, &ksploc));
817 360 : LibmeshPetscCall(KSPSetOperators(ksploc, _hc_sys_h_mat, _hc_sys_h_mat));
818 360 : LibmeshPetscCall(KSPGetPC(ksploc, &pc));
819 360 : LibmeshPetscCall(PCSetType(pc, PCJACOBI));
820 360 : LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit));
821 360 : LibmeshPetscCall(KSPSetOptionsPrefix(ksploc, "h_sys_"));
822 360 : LibmeshPetscCall(KSPSetFromOptions(ksploc));
823 360 : LibmeshPetscCall(KSPSolve(ksploc, _hc_sys_h_rhs, sol));
824 : PetscScalar * xx;
825 360 : LibmeshPetscCall(VecGetArray(sol, &xx));
826 11496 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
827 : {
828 11136 : auto iz_ind = iz - first_node;
829 257160 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
830 : {
831 246024 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
832 246024 : auto h_out = xx[iz_ind * _n_channels + i_ch];
833 246024 : if (h_out < 0)
834 : {
835 0 : mooseError(name(),
836 : " : Calculation of negative Enthalpy h_out = : ",
837 : h_out,
838 : " Axial Level= : ",
839 : iz);
840 : }
841 246024 : _h_soln->set(node_out, h_out);
842 : }
843 : }
844 360 : LibmeshPetscCall(KSPDestroy(&ksploc));
845 360 : LibmeshPetscCall(VecDestroy(&sol));
846 : }
847 : }
848 722 : }
|
