21 params.
addClassDescription(
"Solver class for subchannels in a triangular lattice assembly and " 22 "bare/wire-wrapped fuel pins");
44 PetscErrorCode ierr =
cleanUp();
69 Real standard_area, wire_area, additional_area, wetted_perimeter, displaced_area;
80 auto theta = std::acos(wire_lead_length /
81 std::sqrt(
std::pow(wire_lead_length, 2) +
83 for (
unsigned int iz = 0; iz <
_n_cells + 1; iz++)
85 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
91 Real rod_perimeter = 0.0;
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);
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);
112 additional_area = 0.0;
113 displaced_area = 0.0;
115 wetted_perimeter = rod_perimeter + 0.5 *
libMesh::pi * wire_diameter / std::cos(theta);
119 standard_area =
pitch * (pin_diameter / 2.0 + gap);
120 additional_area = 0.0;
121 displaced_area = (*_displacement_soln)(node)*
pitch;
124 rod_perimeter + 0.5 *
libMesh::pi * wire_diameter / std::cos(theta) +
pitch;
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);
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));
140 auto subchannel_area =
141 standard_area + additional_area + displaced_area - rod_area - wire_area;
144 auto overlapping_pin_area = 0.0;
145 auto overlapping_wetted_perimeter = 0.0;
149 auto pin_1 = gap_pins.first;
150 auto pin_2 = gap_pins.second;
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;
159 if (
pitch < (Radius1 + Radius2))
161 mooseWarning(
" The gap of index : '", i_gap,
" at axial cell ", iz,
" ' is blocked.");
163 (
pitch *
pitch + Radius1 * Radius1 - Radius2 * Radius2) / (2 *
pitch * Radius1);
165 (
pitch *
pitch + Radius2 * Radius2 - Radius1 * Radius1) / (2 *
pitch * Radius2);
166 auto angle1 = 2.0 * acos(cos1);
167 auto angle2 = 2.0 * acos(cos2);
169 overlapping_wetted_perimeter += 0.5 * angle1 * Radius1 + 0.5 * angle2 * Radius2;
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));
177 subchannel_area += overlapping_pin_area;
178 wetted_perimeter += -overlapping_wetted_perimeter;
182 for (
const auto & i_blockage : index_blockage)
184 if (i_ch == i_blockage && (
Z >= z_blockage.front() &&
Z <= z_blockage.back()))
186 subchannel_area *= reduction_blockage[index];
195 for (
unsigned int iz = 0; iz <
_n_cells + 1; iz++)
197 for (
unsigned int i_gap = 0; i_gap <
_n_gaps; i_gap++)
200 auto pin_1 = gap_pins.first;
201 auto pin_2 = gap_pins.second;
207 auto displacement = 0.0;
215 displacement += (*_displacement_soln)(node);
219 displacement = displacement / counter;
221 0.5 * (flat_to_flat - (n_rings - 1) *
pitch * std::sqrt(3.0) -
222 (*_Dpin_soln)(pin_node_1)) +
228 pitch - (*_Dpin_soln)(pin_node_1) / 2.0 - (*
_Dpin_soln)(pin_node_2) / 2.0;
237 for (
unsigned int iz = 1; iz <
_n_cells + 1; iz++)
239 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
249 _aux->solution().close();
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;
262 Real aT, b1T, b2T, cT;
267 auto p_over_d =
pitch / pin_diameter;
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) +
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;
314 cL = aL + b1L * (p_over_d - 1) + b2L *
std::pow((p_over_d - 1), 2.0);
316 cT = aT + b1T * (p_over_d - 1) + b2T *
std::pow((p_over_d - 1), 2.0);
339 cL = aL + b1L * (w_over_d - 1) + b2L *
std::pow((w_over_d - 1), 2.0);
341 cT = aT + b1T * (w_over_d - 1) + b2T *
std::pow((w_over_d - 1), 2.0);
364 cL = aL + b1L * (w_over_d - 1) + b2L *
std::pow((w_over_d - 1), 2.0);
366 cT = aT + b1T * (w_over_d - 1) + b2T *
std::pow((w_over_d - 1), 2.0);
372 if ((wire_diameter != 0.0) && (wire_lead_length != 0.0))
379 ar =
libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
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);
387 cL *= (pw_p / w_perim);
388 cL += wd_l * (3.0 * ar / a_p) * (Dh_i / wire_lead_length) * (Dh_i / wire_diameter);
393 ar =
libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 4.0;
397 cT *=
std::pow((1 + ws_t * (ar / a_p) *
std::pow(std::tan(theta), 2.0)), 1.41);
399 cL *= (1 + ws_l * (ar / a_p) *
std::pow(std::tan(theta), 2.0));
404 ar =
libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
408 cT *=
std::pow((1 + ws_t * (ar / a_p) *
std::pow(std::tan(theta), 2.0)), 1.41);
410 cL *= (1 + ws_l * (ar / a_p) *
std::pow(std::tan(theta), 2.0));
446 unsigned int i_ch = chans.first;
447 unsigned int j_ch = chans.second;
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 +
465 auto avg_hD = 4.0 * (Si + Sj) / (w_perim_i + w_perim_j);
467 0.5 * (((*_mdot_soln)(node_in_i) + (*
_mdot_soln)(node_in_j)) / (Si_in + Sj_in) +
469 auto Re = avg_massflux * avg_hD / avg_mu;
474 auto ReT = 10000.0 *
std::pow(10.0, 0.7 * (
pitch / pin_diameter - 1));
478 (wire_lead_length != 0) && (wire_diameter != 0))
481 auto theta = std::acos(wire_lead_length /
482 std::sqrt(
std::pow(wire_lead_length, 2) +
484 auto Ar1 =
libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
491 Cm = 0.077 *
std::pow((
pitch - pin_diameter) / pin_diameter, -0.5);
495 Cm = 0.14 *
std::pow((
pitch - pin_diameter) / pin_diameter, -0.5);
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)) *
507 beta = Cm *
std::pow(Ar1 /
A1, 0.5) * std::tan(theta);
510 else if ((wire_lead_length == 0) && (wire_diameter == 0))
517 auto k = (1 / S_total) *
522 auto cp = (1 / S_total) *
527 auto Pr = avg_mu *
cp /
k;
528 auto Pr_t = Pr * (Re / gamma) * std::sqrt(
f / 8.0);
530 auto L_x = sf *
delta;
531 auto lamda = gap / L_x;
532 auto a_x = 1.0 - 2.0 * lamda * lamda /
libMesh::pi;
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);
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;
540 beta = freq_factor * (rod_mixing + axial_mixing) *
std::pow(Re, -
b / 2.0);
542 mooseAssert(beta >= 0,
543 "beta should be positive for the inner gaps, or zero for the edge gaps, because this " 545 "explicitly in the computeh method.");
580 double heat_rate_in = 0.0;
581 double heat_rate_out = 0.0;
586 heat_rate_out += factor * (*_q_prime_soln)(node_out);
587 heat_rate_in += factor * (*_q_prime_soln)(node_in);
589 return (heat_rate_in + heat_rate_out) * dz / 2.0;
605 unsigned int last_node = (iblock + 1) *
_block_size;
606 unsigned int first_node = iblock *
_block_size + 1;
614 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
621 name(),
" : Calculation of negative Enthalpy h_out = : ", h_out,
" Axial Level= : ", 0);
629 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
632 auto dz = z_grid[iz] - z_grid[iz - 1];
633 Real gedge_ave = 0.0;
636 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
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;
648 gedge_ave = mdot_sum / si_sum;
650 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
654 auto mdot_in = (*_mdot_soln)(node_in);
655 auto h_in = (*_h_soln)(node_in);
656 auto volume = dz * (*_S_flow_soln)(node_in);
657 auto mdot_out = (*_mdot_soln)(node_out);
660 Real sumWijPrimeDhij = 0.0;
668 Real sweep_enthalpy = 0.0;
671 (wire_diameter != 0.0) && (wire_lead_length != 0.0))
678 auto w = pin_diameter + gap;
680 std::acos(wire_lead_length /
681 std::sqrt(
std::pow(wire_lead_length, 2) +
684 auto Si = (*_S_flow_soln)(node_in);
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);
696 auto hD = 4.0 * Si / w_perim;
697 auto Re = massflux * hD /
mu;
699 auto Ar2 =
libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 4.0;
706 Cs = 0.033 *
std::pow(wire_lead_length / pin_diameter, 0.3);
710 Cs = 0.75 *
std::pow(wire_lead_length / pin_diameter, 0.3);
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)) *
722 auto beta = Cs *
std::pow(Ar2 /
A2, 0.5) * std::tan(theta);
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);
732 unsigned int counter = 0;
736 unsigned int ii_ch = chans.first;
738 unsigned int jj_ch = chans.second;
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);
749 sumWijPrimeDhij +=
_WijPrime(i_gap, iz) * (2 * (*_h_soln)(node_in) -
769 dist_ij =
pitch / std::sqrt(3);
776 0.66 * (
pitch / pin_diameter) *
780 e_cond += 0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
785 e_cond += -0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
795 auto Si = (*_S_flow_soln)(node_in_i);
796 auto dist_ij = z_grid[iz] - z_grid[iz - 1];
798 e_cond += 0.5 * (thcon_i + thcon_j) * Si * ((*
_T_soln)(node_in_j) - (*
_T_soln)(node_in_i)) /
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 *
817 (mdot_in * h_in - sumWijh - sumWijPrimeDhij + added_enthalpy + e_cond + sweep_enthalpy +
823 " : Calculation of negative Enthalpy h_out = : ",
852 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
857 auto iz_ind = iz - first_node;
859 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
863 auto S_in = (*_S_flow_soln)(node_in);
864 auto S_out = (*_S_flow_soln)(node_out);
866 auto volume = dz * S_interp;
868 PetscScalar Pe = 0.5;
872 auto w_perim_in = (*_w_perim_soln)(node_in);
873 auto w_perim_out = (*_w_perim_soln)(node_out);
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));
890 if (iz == first_node)
892 PetscScalar value_vec_tt =
901 PetscInt col_tt = i_ch +
_n_channels * (iz_ind - 1);
903 LibmeshPetscCall(MatSetValues(
911 LibmeshPetscCall(MatSetValues(
915 PetscScalar rho_old_interp =
917 PetscScalar h_old_interp =
919 PetscScalar value_vec_tt =
_TR * rho_old_interp * h_old_interp *
volume /
_dt;
925 if (iz == first_node)
928 PetscScalar value_at =
alpha * (*_mdot_soln)(node_in) * (*
_h_soln)(node_in);
932 value_at =
alpha * (*_mdot_soln)(node_out) - (1 -
alpha) * (*_mdot_soln)(node_in);
934 LibmeshPetscCall(MatSetValues(
939 LibmeshPetscCall(MatSetValues(
942 else if (iz == last_node)
945 PetscScalar value_at = 1.0 * (*_mdot_soln)(node_out);
947 LibmeshPetscCall(MatSetValues(
950 value_at = -1.0 * (*_mdot_soln)(node_in);
952 LibmeshPetscCall(MatSetValues(
960 PetscScalar value_at = -
alpha * (*_mdot_soln)(node_in);
962 LibmeshPetscCall(MatSetValues(
965 value_at =
alpha * (*_mdot_soln)(node_out) - (1 -
alpha) * (*_mdot_soln)(node_in);
967 LibmeshPetscCall(MatSetValues(
972 LibmeshPetscCall(MatSetValues(
981 auto diff_center = K_center / (cp_center + 1e-15);
983 if (iz == first_node)
993 auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
994 auto diff_top = K_top / (cp_top + 1e-15);
1008 PetscScalar value_at = diff_up * S_up / dz_up + diff_down * S_down / dz_down;
1009 LibmeshPetscCall(MatSetValues(
1013 value_at = 1.0 * diff_down * S_down / dz_down * (*_h_soln)(node_bottom);
1019 value_at = -diff_up * S_up / dz_up;
1020 LibmeshPetscCall(MatSetValues(
1023 else if (iz == last_node)
1030 auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
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);
1039 PetscScalar value_at = diff_down * S_down / dz_down;
1040 LibmeshPetscCall(MatSetValues(
1045 value_at = -diff_down * S_down / dz_down;
1046 LibmeshPetscCall(MatSetValues(
1064 auto diff_bottom = K_bottom / (cp_bottom + 1e-15);
1065 auto diff_top = K_top / (cp_top + 1e-15);
1079 PetscScalar value_at = diff_up * S_up / dz_up + diff_down * S_down / dz_down;
1080 LibmeshPetscCall(MatSetValues(
1085 value_at = -diff_down * S_down / dz_down;
1086 LibmeshPetscCall(MatSetValues(
1091 value_at = -diff_up * S_up / dz_up;
1092 LibmeshPetscCall(MatSetValues(
1097 unsigned int counter = 0;
1098 unsigned int cross_index = iz;
1103 unsigned int ii_ch = chans.first;
1104 unsigned int jj_ch = chans.second;
1109 if (
_Wij(i_gap, cross_index) > 0.0)
1111 if (iz == first_node)
1113 h_star = (*_h_soln)(node_in_i);
1114 PetscScalar value_vec_ct = -1.0 *
alpha *
1116 _Wij(i_gap, cross_index) * h_star;
1117 PetscInt row_vec_ct = i_ch +
_n_channels * iz_ind;
1118 LibmeshPetscCall(VecSetValues(
1124 _Wij(i_gap, cross_index);
1126 PetscInt col_ct = ii_ch +
_n_channels * (iz_ind - 1);
1127 LibmeshPetscCall(MatSetValues(
1130 PetscScalar value_ct = (1.0 -
alpha) *
1132 _Wij(i_gap, cross_index);
1135 LibmeshPetscCall(MatSetValues(
1138 else if (
_Wij(i_gap, cross_index) < 0.0)
1140 if (iz == first_node)
1142 h_star = (*_h_soln)(node_in_j);
1143 PetscScalar value_vec_ct = -1.0 *
alpha *
1145 _Wij(i_gap, cross_index) * h_star;
1146 PetscInt row_vec_ct = i_ch +
_n_channels * iz_ind;
1147 LibmeshPetscCall(VecSetValues(
1153 _Wij(i_gap, cross_index);
1155 PetscInt col_ct = jj_ch +
_n_channels * (iz_ind - 1);
1156 LibmeshPetscCall(MatSetValues(
1159 PetscScalar value_ct = (1.0 -
alpha) *
1161 _Wij(i_gap, cross_index);
1164 LibmeshPetscCall(MatSetValues(
1169 if (iz == first_node)
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;
1181 PetscScalar value_center_ct = 2.0 *
alpha *
_WijPrime(i_gap, cross_index);
1183 PetscInt col_ct = i_ch +
_n_channels * (iz_ind - 1);
1184 LibmeshPetscCall(MatSetValues(
1187 PetscScalar value_left_ct = -1.0 *
alpha *
_WijPrime(i_gap, cross_index);
1190 LibmeshPetscCall(MatSetValues(
1193 PetscScalar value_right_ct = -1.0 *
alpha *
_WijPrime(i_gap, cross_index);
1196 LibmeshPetscCall(MatSetValues(
1199 PetscScalar value_center_ct = 2.0 * (1.0 -
alpha) *
_WijPrime(i_gap, cross_index);
1202 LibmeshPetscCall(MatSetValues(
1205 PetscScalar value_left_ct = -1.0 * (1.0 -
alpha) *
_WijPrime(i_gap, cross_index);
1208 LibmeshPetscCall(MatSetValues(
1211 PetscScalar value_right_ct = -1.0 * (1.0 -
alpha) *
_WijPrime(i_gap, cross_index);
1214 LibmeshPetscCall(MatSetValues(
1233 dist_ij =
pitch / std::sqrt(3);
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);
1245 0.66 * (
pitch / pin_diameter) *
1248 auto base_value = harm_A * shape_factor * Sij / dist_ij;
1249 auto neg_base_value = -1.0 * base_value;
1253 LibmeshPetscCall(MatSetValues(
1258 LibmeshPetscCall(MatSetValues(
1263 LibmeshPetscCall(MatSetValues(
1268 LibmeshPetscCall(MatSetValues(
1274 Real gedge_ave = 0.0;
1275 Real mdot_sum = 0.0;
1277 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
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;
1289 gedge_ave = mdot_sum / si_sum;
1291 PetscScalar sweep_enthalpy = 0.0;
1293 (wire_diameter != 0.0) && (wire_lead_length != 0.0))
1300 auto w = pin_diameter + gap;
1302 std::acos(wire_lead_length /
1303 std::sqrt(
std::pow(wire_lead_length, 2) +
1305 auto Sij = dz * gap;
1306 auto Si = (*_S_flow_soln)(node_in);
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);
1318 auto hD = 4.0 * Si / w_perim;
1319 auto Re = massflux * hD /
mu;
1321 auto Ar2 =
libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 4.0;
1328 Cs = 0.033 *
std::pow(wire_lead_length / pin_diameter, 0.3);
1332 Cs = 0.75 *
std::pow(wire_lead_length / pin_diameter, 0.3);
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)) *
1344 auto beta = Cs *
std::pow(Ar2 /
A2, 0.5) * std::tan(theta);
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);
1352 if (iz == first_node)
1355 PetscScalar value_hs = -sweep_enthalpy;
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(
1365 PetscInt col_sh_l = sweep_in +
_n_channels * (iz_ind - 1);
1366 PetscScalar neg_sweep_in = -1.0 * wsweep_in;
1367 LibmeshPetscCall(MatSetValues(
1375 PetscInt row_vec_ht = i_ch +
_n_channels * iz_ind;
1393 LibmeshPetscCall(MatAssemblyBegin(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1394 LibmeshPetscCall(MatAssemblyEnd(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1398 LibmeshPetscCall(MatAssemblyBegin(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1399 LibmeshPetscCall(MatAssemblyEnd(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1402 LibmeshPetscCall(MatAssemblyBegin(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1403 LibmeshPetscCall(MatAssemblyEnd(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1406 LibmeshPetscCall(MatAssemblyBegin(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1407 LibmeshPetscCall(MatAssemblyEnd(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1410 LibmeshPetscCall(MatAssemblyBegin(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1411 LibmeshPetscCall(MatAssemblyEnd(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1414 LibmeshPetscCall(MatAssemblyBegin(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1415 LibmeshPetscCall(MatAssemblyEnd(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1418 LibmeshPetscCall(MatAssemblyBegin(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1419 LibmeshPetscCall(MatAssemblyEnd(
_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
1421 _console <<
"Block: " << iblock <<
" - Enthalpy conservation matrix assembled" << std::endl;
1438 LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
1440 LibmeshPetscCall(KSPGetPC(ksploc, &pc));
1441 LibmeshPetscCall(PCSetType(pc, PCJACOBI));
1443 LibmeshPetscCall(KSPSetOptionsPrefix(ksploc,
"h_sys_"));
1444 LibmeshPetscCall(KSPSetFromOptions(ksploc));
1448 LibmeshPetscCall(VecGetArray(sol, &xx));
1449 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
1451 auto iz_ind = iz - first_node;
1452 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1459 " : Calculation of negative Enthalpy h_out = : ",
1464 _h_soln->set(node_out, h_out);
1467 LibmeshPetscCall(KSPDestroy(&ksploc));
1468 LibmeshPetscCall(VecDestroy(&sol));
Vec _hc_radial_heat_conduction_rhs
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.
Mat _hc_radial_heat_conduction_mat
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
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.
SubChannelMesh & _subchannel_mesh
virtual Real computeFrictionFactor(FrictionStruct friction_args) override
Computes the axial friction factor for the sodium coolant and for each subchannel.
TriSubChannelMesh & _tri_sch_mesh
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
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.
Vec _hc_time_derivative_rhs
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$.
static InputParameters validParams()
Mat _hc_sweep_enthalpy_mat
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.
Vec _hc_cross_derivative_rhs
void mooseWarning(Args &&... args) const
Real computeAddedHeatDuct(unsigned int i_ch, unsigned int iz)
Function that computes the heat flux added by the duct.
Mat _hc_axial_heat_conduction_mat
Vec _hc_sweep_enthalpy_rhs
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
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
Real f(Real x)
Test function for Brents method.
std::unique_ptr< SolutionHandle > _q_prime_soln
static const std::string mu
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.
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
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
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.
Vec _hc_advective_derivative_rhs
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.
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)
Vec _hc_axial_heat_conduction_rhs
TriSubChannel1PhaseProblem(const InputParameters ¶ms)
MooseUnits pow(const MooseUnits &, int)
virtual ~TriSubChannel1PhaseProblem()
std::unique_ptr< SolutionHandle > _P_soln
static const std::string k
std::unique_ptr< SolutionHandle > _w_perim_soln
static const std::string cL
virtual Real computeBeta(unsigned int i_gap, unsigned int iz) override
Computes turbulent mixing coefficient.
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.
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.