Skip to content

Commit 1046d4b

Browse files
kyriv1980kyriv1980
committed
address reviewer comments #30049 
1 parent 855f30b commit 1046d4b

File tree

10 files changed

+282
-244
lines changed

10 files changed

+282
-244
lines changed

‎modules/subchannel/doc/content/source/scmclosures/SCMMixingChengTodreas.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66

77
!! Intentional comment to provide extra spacing
88

9-
This closure class is used to model the turbulent mixing coefficient $\beta$ using the Cheng and Todreas corelations. Specifically this closure model applies to triangular assemblies with wire-wrapped pins. Citation:
9+
This closure class is used to model the turbulent mixing coefficient $\beta$ using the Cheng and Todreas correlations. Specifically this closure model applies to triangular assemblies with wire-wrapped pins. The implementation was based on:
1010

1111
- Hydrodynamic models and correlations for bare and wire-wrapped hexagonal rod bundles—bundle friction factors, subchannel friction factors and mixing parameters, Cheng and Todreas [!cite](cheng1986hydrodynamic).
1212

‎modules/subchannel/doc/content/source/scmclosures/SCMMixingConstantBeta.md‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,10 +18,12 @@ It is important to note that the mixing coefficient is simply a tuning parameter
1818

1919
After calibrating the turbulent diffusion coefficient $\beta$ we turned our attention to the turbulent modeling parameter $C_{T}$. This is a tuning parameter that informs on how much momentum is transferred/diffused between subchannels, due to turbulence. The CNEN 4x4 test [!cite](Marinelli) performed at Studsvik laboratory for studying the flow mixing effect between adjacent subchannels was chosen to tune this parameter. This experiment consists in velocity and temperature measurements taken at the outlet of a 16-pin assembly test section. Analysis of the velocity distribution at the exit of the assembly can be used to calibrate the turbulent parameter $C_{T}$.
2020

21-
For quadrilateral assemblies: $C_{T} = 2.6$, $\beta = 0.006$ [!cite](kyriakopoulos2022development).
21+
For quadrilateral assemblies, the calibration values computed were: $C_{T} = 2.6$, $\beta = 0.006$ [!cite](kyriakopoulos2022development).
2222

2323
!syntax parameters /SCMClosures/SCMMixingConstantBeta
2424

2525
!syntax inputs /SCMClosures/SCMMixingConstantBeta
2626

2727
!syntax children /SCMClosures/SCMMixingConstantBeta
28+
29+
!bibtex bibliography

‎modules/subchannel/doc/content/source/scmclosures/SCMMixingKimAndChung.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66

77
!! Intentional comment to provide extra spacing
88

9-
This closure class is used to model the turbulent mixing coefficient $\beta$ using the Kim and Chung corelations. Specifically this closure model applies to triangular and quadrilateral assemblies with bare pins. Citations:
9+
This closure class is used to model the turbulent mixing coefficient $\beta$ using the Kim and Chung correlations. Specifically this closure model applies to triangular and quadrilateral assemblies with bare pins. The implementation followed:
1010

1111
- A scale analysis of the turbulent mixing rate for various Prandtl number flow fields in rod bundles eq 25,Kim and Chung (2001) [!cite](kim2001scale).
1212
- Modeling of flow blockage in a liquid metal-cooled reactor subassembly with a subchannel analysis code eq 19, Jeong et. al (2005)[!cite](jeong2005modeling).

‎modules/subchannel/include/scmclosures/SCMMixingChengTodreas.h‎

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -24,12 +24,17 @@ class SCMMixingChengTodreas : public SCMMixingClosureBase
2424

2525
SCMMixingChengTodreas(const InputParameters & parameters);
2626

27-
virtual Real computeMixingParameter(const unsigned int & i_gap,
28-
const unsigned int & iz,
29-
const bool & sweep_flow) const override;
27+
Real computeMixingParameter(const unsigned int i_gap,
28+
const unsigned int iz,
29+
const bool sweep_flow) const override;
3030

3131
/// Keep track of the lattice type
3232
bool _is_tri_lattice;
3333
/// Pointer to the tri lattice mesh
3434
const TriSubChannelMesh * const _tri_sch_mesh;
35+
36+
SolutionHandle _S_soln;
37+
SolutionHandle _mdot_soln;
38+
SolutionHandle _w_perim_soln;
39+
SolutionHandle _mu_soln;
3540
};

‎modules/subchannel/include/scmclosures/SCMMixingClosureBase.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ class SCMMixingClosureBase : public SCMClosureBase
2727
/// @brief Computes the turbulent mixing coefficient for the local conditions around gap(i_gap) and axial level(iz)
2828
/// @param i_gap and @param iz
2929
/// @return the mixing coefficient (beta)
30-
virtual Real computeMixingParameter(const unsigned int & i_gap,
31-
const unsigned int & iz,
32-
const bool & sweep_flow = false) const = 0;
30+
virtual Real computeMixingParameter(const unsigned int i_gap,
31+
const unsigned int iz,
32+
const bool sweep_flow = false) const = 0;
3333
};

‎modules/subchannel/include/scmclosures/SCMMixingConstantBeta.h‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -23,11 +23,11 @@ class SCMMixingConstantBeta : public SCMMixingClosureBase
2323

2424
SCMMixingConstantBeta(const InputParameters & parameters);
2525

26-
virtual Real computeMixingParameter(const unsigned int & i_gap,
27-
const unsigned int & iz,
28-
const bool & sweep_flow) const override;
26+
virtual Real computeMixingParameter(const unsigned int i_gap,
27+
const unsigned int iz,
28+
const bool sweep_flow) const override;
2929

3030
protected:
3131
/// Turbulent mixing parameter
32-
const Real & _beta;
32+
const Real _beta;
3333
};

‎modules/subchannel/include/scmclosures/SCMMixingKimAndChung.h‎

Lines changed: 28 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -10,12 +10,12 @@
1010
#pragma once
1111

1212
#include "SCMMixingClosureBase.h"
13-
#include "TriSubChannelMesh.h"
14-
#include "QuadSubChannelMesh.h"
13+
14+
class SubChannelMesh;
1515

1616
/**
1717
* Class that calculates the turbulent mixing coefficient based on the Kim and Chung (2001)
18-
* corellation (eq 25) It is used for both quad and tri lattices with bare fuel pins.
18+
* correllation (eq 25). It is used for both quad and tri lattices with bare fuel pins.
1919
*/
2020
class SCMMixingKimAndChung : public SCMMixingClosureBase
2121
{
@@ -24,22 +24,34 @@ class SCMMixingKimAndChung : public SCMMixingClosureBase
2424

2525
SCMMixingKimAndChung(const InputParameters & parameters);
2626

27-
virtual Real computeMixingParameter(const unsigned int & i_gap,
28-
const unsigned int & iz,
29-
const bool & sweep_flow) const override;
27+
virtual Real computeMixingParameter(const unsigned int i_gap,
28+
const unsigned int iz,
29+
const bool sweep_flow) const override;
3030

3131
protected:
32-
Real computeTriLatticeMixingParameter(const unsigned int & i_gap,
33-
const unsigned int & iz,
34-
const bool & sweep_flow) const;
35-
Real computeQuadLatticeMixingParameter(const unsigned int & i_gap,
36-
const unsigned int & iz,
37-
const bool & sweep_flow) const;
32+
Real computeTriLatticeMixingParameter(const unsigned int i_gap,
33+
const unsigned int iz,
34+
const bool sweep_flow) const;
35+
36+
Real computeQuadLatticeMixingParameter(const unsigned int i_gap,
37+
const unsigned int iz,
38+
const bool sweep_flow) const;
39+
40+
Real computeLatticeMixingParameter(const unsigned int i_gap,
41+
const unsigned int iz,
42+
const Real delta,
43+
const Real sf,
44+
const bool sweep_flow) const;
3845

3946
/// Keep track of the lattice type
4047
bool _is_tri_lattice;
41-
/// Pointer to the tri lattice mesh
42-
const TriSubChannelMesh * const _tri_sch_mesh;
43-
/// Pointer to the quad lattice mesh
44-
const QuadSubChannelMesh * const _quad_sch_mesh;
48+
/// Pointer to the base subchannel mesh
49+
const SubChannelMesh & _sch_mesh;
50+
51+
SolutionHandle _S_soln;
52+
SolutionHandle _mdot_soln;
53+
SolutionHandle _w_perim_soln;
54+
SolutionHandle _mu_soln;
55+
SolutionHandle _P_soln;
56+
SolutionHandle _T_soln;
4557
};

‎modules/subchannel/src/scmclosures/SCMMixingChengTodreas.C‎

Lines changed: 99 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -23,57 +23,75 @@ SCMMixingChengTodreas::validParams()
2323
SCMMixingChengTodreas::SCMMixingChengTodreas(const InputParameters & parameters)
2424
: SCMMixingClosureBase(parameters),
2525
_is_tri_lattice(dynamic_cast<const TriSubChannelMesh *>(&_subchannel_mesh) != nullptr),
26-
_tri_sch_mesh(dynamic_cast<const TriSubChannelMesh *>(&_subchannel_mesh))
26+
_tri_sch_mesh(dynamic_cast<const TriSubChannelMesh *>(&_subchannel_mesh)),
27+
_S_soln(_subproblem.getVariable(0, "S")),
28+
_mdot_soln(_subproblem.getVariable(0, "mdot")),
29+
_w_perim_soln(_subproblem.getVariable(0, "w_perim")),
30+
_mu_soln(_subproblem.getVariable(0, "mu"))
2731
{
2832
}
2933

3034
Real
31-
SCMMixingChengTodreas::computeMixingParameter(const unsigned int & i_gap,
32-
const unsigned int & iz,
33-
const bool & sweep_flow) const
35+
SCMMixingChengTodreas::computeMixingParameter(const unsigned int i_gap,
36+
const unsigned int iz,
37+
const bool sweep_flow) const
3438
{
3539
if (!_is_tri_lattice)
3640
mooseError("This corelation applies only for triangular assemblies");
37-
auto beta = std::numeric_limits<double>::quiet_NaN();
38-
const Real & pitch = _subchannel_mesh.getPitch();
39-
const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
40-
const Real & wire_lead_length = _tri_sch_mesh->getWireLeadLength();
41-
const Real & wire_diameter = _tri_sch_mesh->getWireDiameter();
41+
42+
Real beta = std::numeric_limits<double>::quiet_NaN();
43+
44+
const Real pitch = _subchannel_mesh.getPitch();
45+
const Real pin_diameter = _subchannel_mesh.getPinDiameter();
46+
47+
const Real wire_lead_length = _tri_sch_mesh->getWireLeadLength();
48+
const Real wire_diameter = _tri_sch_mesh->getWireDiameter();
49+
4250
if (wire_lead_length == 0 && wire_diameter == 0)
4351
mooseError("This corelation applies only for wire-wrapped assemblies");
44-
auto S_soln = SolutionHandle(_subproblem.getVariable(0, "S"));
45-
auto mdot_soln = SolutionHandle(_subproblem.getVariable(0, "mdot"));
46-
auto w_perim_soln = SolutionHandle(_subproblem.getVariable(0, "w_perim"));
47-
auto mu_soln = SolutionHandle(_subproblem.getVariable(0, "mu"));
48-
auto chans = _subchannel_mesh.getGapChannels(i_gap);
49-
auto Nr = _tri_sch_mesh->getNumOfRings();
50-
unsigned int i_ch = chans.first;
51-
unsigned int j_ch = chans.second;
52-
auto subch_type_i = _subchannel_mesh.getSubchannelType(i_ch);
53-
auto subch_type_j = _subchannel_mesh.getSubchannelType(j_ch);
54-
auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
55-
auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
56-
auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
57-
auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
58-
auto Si_in = (S_soln)(node_in_i);
59-
auto Sj_in = (S_soln)(node_in_j);
60-
auto Si_out = (S_soln)(node_out_i);
61-
auto Sj_out = (S_soln)(node_out_j);
62-
auto S_total = Si_in + Sj_in + Si_out + Sj_out;
63-
auto Si = 0.5 * (Si_in + Si_out);
64-
auto Sj = 0.5 * (Sj_in + Sj_out);
65-
auto w_perim_i = 0.5 * ((w_perim_soln)(node_in_i) + (w_perim_soln)(node_out_i));
66-
auto w_perim_j = 0.5 * ((w_perim_soln)(node_in_j) + (w_perim_soln)(node_out_j));
67-
auto avg_mu = (1 / S_total) * ((mu_soln)(node_out_i)*Si_out + (mu_soln)(node_in_i)*Si_in +
68-
(mu_soln)(node_out_j)*Sj_out + (mu_soln)(node_in_j)*Sj_in);
69-
auto avg_hD = 4.0 * (Si + Sj) / (w_perim_i + w_perim_j);
70-
auto avg_massflux =
71-
0.5 * (((mdot_soln)(node_in_i) + (mdot_soln)(node_in_j)) / (Si_in + Sj_in) +
72-
((mdot_soln)(node_out_i) + (mdot_soln)(node_out_j)) / (Si_out + Sj_out));
73-
auto Re = avg_massflux * avg_hD / avg_mu;
52+
53+
const auto chans = _subchannel_mesh.getGapChannels(i_gap);
54+
const unsigned int i_ch = chans.first;
55+
const unsigned int j_ch = chans.second;
56+
57+
const unsigned int Nr = _tri_sch_mesh->getNumOfRings();
58+
59+
const auto subch_type_i = _subchannel_mesh.getSubchannelType(i_ch);
60+
const auto subch_type_j = _subchannel_mesh.getSubchannelType(j_ch);
61+
62+
const Node * const node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
63+
const Node * const node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
64+
const Node * const node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
65+
const Node * const node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
66+
67+
const Real Si_in = _S_soln(node_in_i);
68+
const Real Sj_in = _S_soln(node_in_j);
69+
const Real Si_out = _S_soln(node_out_i);
70+
const Real Sj_out = _S_soln(node_out_j);
71+
72+
const Real S_total = Si_in + Sj_in + Si_out + Sj_out;
73+
const Real Si = 0.5 * (Si_in + Si_out);
74+
const Real Sj = 0.5 * (Sj_in + Sj_out);
75+
76+
const Real w_perim_i = 0.5 * (_w_perim_soln(node_in_i) + _w_perim_soln(node_out_i));
77+
const Real w_perim_j = 0.5 * (_w_perim_soln(node_in_j) + _w_perim_soln(node_out_j));
78+
79+
const Real avg_mu =
80+
(1.0 / S_total) * (_mu_soln(node_out_i) * Si_out + _mu_soln(node_in_i) * Si_in +
81+
_mu_soln(node_out_j) * Sj_out + _mu_soln(node_in_j) * Sj_in);
82+
83+
const Real avg_hD = 4.0 * (Si + Sj) / (w_perim_i + w_perim_j);
84+
85+
const Real avg_massflux =
86+
0.5 * ((_mdot_soln(node_in_i) + _mdot_soln(node_in_j)) / (Si_in + Sj_in) +
87+
(_mdot_soln(node_out_i) + _mdot_soln(node_out_j)) / (Si_out + Sj_out));
88+
89+
const Real Re = avg_massflux * avg_hD / avg_mu;
90+
7491
// Calculation of flow regime
75-
auto ReL = 320.0 * std::pow(10.0, pitch / pin_diameter - 1);
76-
auto ReT = 10000.0 * std::pow(10.0, 0.7 * (pitch / pin_diameter - 1));
92+
const Real ReL = 320.0 * std::pow(10.0, pitch / pin_diameter - 1.0);
93+
const Real ReT = 10000.0 * std::pow(10.0, 0.7 * (pitch / pin_diameter - 1.0));
94+
7795
// Calculation of Turbulent Crossflow for wire-wrapped triangular assemblies. Cheng &
7896
// Todreas (1986).
7997
// INNER SUBCHANNELS
@@ -82,21 +100,25 @@ SCMMixingChengTodreas::computeMixingParameter(const unsigned int & i_gap,
82100
{
83101
// Calculation of geometric parameters
84102
// wire angle
85-
auto theta =
103+
const Real theta =
86104
std::acos(wire_lead_length /
87105
std::sqrt(Utility::pow<2>(wire_lead_length) +
88106
Utility::pow<2>(libMesh::pi * (pin_diameter + wire_diameter))));
107+
89108
// projected area of wire on subchannel
90-
auto Ar1 = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
109+
const Real Ar1 = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
110+
91111
// bare subchannel flow area
92-
auto A1prime = (std::sqrt(3.0) / 4.0) * Utility::pow<2>(pitch) -
93-
libMesh::pi * Utility::pow<2>(pin_diameter) / 8.0;
112+
const Real A1prime = (std::sqrt(3.0) / 4.0) * Utility::pow<2>(pitch) -
113+
libMesh::pi * Utility::pow<2>(pin_diameter) / 8.0;
114+
94115
// wire-wrapped subchannel flow area
95-
auto A1 = A1prime - libMesh::pi * Utility::pow<2>(wire_diameter) / 8.0 / std::cos(theta);
116+
const Real A1 = A1prime - libMesh::pi * Utility::pow<2>(wire_diameter) / 8.0 / std::cos(theta);
117+
96118
// empirical constant for mixing parameter
97-
auto Cm = 0.0;
98-
auto CmL_constant = 0.0;
99-
auto CmT_constant = 0.0;
119+
Real Cm = 0.0;
120+
Real CmL_constant = 0.0;
121+
Real CmT_constant = 0.0;
100122

101123
if (Nr == 1)
102124
{
@@ -109,8 +131,8 @@ SCMMixingChengTodreas::computeMixingParameter(const unsigned int & i_gap,
109131
CmL_constant = 0.077;
110132
}
111133

112-
auto CmT = CmT_constant * std::pow((pitch - pin_diameter) / pin_diameter, -0.5);
113-
auto CmL = CmL_constant * std::pow((pitch - pin_diameter) / pin_diameter, -0.5);
134+
const Real CmT = CmT_constant * std::pow((pitch - pin_diameter) / pin_diameter, -0.5);
135+
const Real CmL = CmL_constant * std::pow((pitch - pin_diameter) / pin_diameter, -0.5);
114136

115137
if (Re < ReL)
116138
{
@@ -122,10 +144,11 @@ SCMMixingChengTodreas::computeMixingParameter(const unsigned int & i_gap,
122144
}
123145
else
124146
{
125-
auto psi = (std::log(Re) - std::log(ReL)) / (std::log(ReT) - std::log(ReL));
126-
auto gamma = 2.0 / 3.0;
147+
const Real psi = (std::log(Re) - std::log(ReL)) / (std::log(ReT) - std::log(ReL));
148+
const Real gamma = 2.0 / 3.0;
127149
Cm = CmL + (CmT - CmL) * std::pow(psi, gamma);
128150
}
151+
129152
// mixing parameter
130153
beta = Cm * std::sqrt(Ar1 / A1) * std::tan(theta);
131154
}
@@ -134,23 +157,30 @@ SCMMixingChengTodreas::computeMixingParameter(const unsigned int & i_gap,
134157
(subch_type_j == EChannelType::CORNER || subch_type_j == EChannelType::EDGE) &&
135158
(wire_lead_length != 0) && (wire_diameter != 0))
136159
{
137-
auto theta =
160+
const Real theta =
138161
std::acos(wire_lead_length /
139162
std::sqrt(Utility::pow<2>(wire_lead_length) +
140163
Utility::pow<2>(libMesh::pi * (pin_diameter + wire_diameter))));
164+
141165
// Calculation of geometric parameters
142166
// distance from pin surface to duct
143-
auto dpgap = _tri_sch_mesh->getDuctToPinGap();
167+
const Real dpgap = _tri_sch_mesh->getDuctToPinGap();
168+
144169
// Edge pitch parameter defined as pin diameter plus distance to duct wall
145-
auto w = pin_diameter + dpgap;
146-
auto Ar2 = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 4.0;
147-
auto A2prime =
170+
const Real w = pin_diameter + dpgap;
171+
172+
const Real Ar2 = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 4.0;
173+
174+
const Real A2prime =
148175
pitch * (w - pin_diameter / 2.0) - libMesh::pi * Utility::pow<2>(pin_diameter) / 8.0;
149-
auto A2 = A2prime - libMesh::pi * Utility::pow<2>(wire_diameter) / 8.0 / std::cos(theta);
176+
177+
const Real A2 = A2prime - libMesh::pi * Utility::pow<2>(wire_diameter) / 8.0 / std::cos(theta);
178+
150179
// empirical constant for mixing parameter
151-
auto Cs = 0.0;
152-
auto CsL_constant = 0.0;
153-
auto CsT_constant = 0.0;
180+
Real Cs = 0.0;
181+
Real CsL_constant = 0.0;
182+
Real CsT_constant = 0.0;
183+
154184
if (Nr == 1)
155185
{
156186
CsT_constant = 0.6;
@@ -161,8 +191,9 @@ SCMMixingChengTodreas::computeMixingParameter(const unsigned int & i_gap,
161191
CsT_constant = 0.75;
162192
CsL_constant = 0.413;
163193
}
164-
auto CsL = CsL_constant * std::pow(wire_lead_length / pin_diameter, 0.3);
165-
auto CsT = CsT_constant * std::pow(wire_lead_length / pin_diameter, 0.3);
194+
195+
const Real CsL = CsL_constant * std::pow(wire_lead_length / pin_diameter, 0.3);
196+
const Real CsT = CsT_constant * std::pow(wire_lead_length / pin_diameter, 0.3);
166197

167198
if (Re < ReL)
168199
{
@@ -174,15 +205,17 @@ SCMMixingChengTodreas::computeMixingParameter(const unsigned int & i_gap,
174205
}
175206
else
176207
{
177-
auto psi = (std::log(Re) - std::log(ReL)) / (std::log(ReT) - std::log(ReL));
178-
auto gamma = 2.0 / 3.0;
208+
const Real psi = (std::log(Re) - std::log(ReL)) / (std::log(ReT) - std::log(ReL));
209+
const Real gamma = 2.0 / 3.0;
179210
Cs = CsL + (CsT - CsL) * std::pow(psi, gamma);
180211
}
212+
181213
// Calculation of turbulent mixing parameter used for sweep flow only in enthalpy calculation
182214
if (sweep_flow)
183215
beta = Cs * std::sqrt(Ar2 / A2) * std::tan(theta);
184216
else
185217
beta = 0.0;
186218
}
219+
187220
return beta;
188221
}

0 commit comments

Comments
 (0)