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 "SCMMixingKimAndChung.h"
      11             : #include "SCMFrictionClosureBase.h"
      12             : #include "SinglePhaseFluidProperties.h"
      13             : #include "TriSubChannelMesh.h"
      14             : #include "QuadSubChannelMesh.h"
      15             : 
      16             : registerMooseObject("SubChannelApp", SCMMixingKimAndChung);
      17             : 
      18             : InputParameters
      19          57 : SCMMixingKimAndChung::validParams()
      20             : {
      21          57 :   InputParameters params = SCMMixingClosureBase::validParams();
      22          57 :   params.addClassDescription(
      23             :       "Class that models the turbulent mixing coefficient using the Kim and Chung correlations.");
      24          57 :   return params;
      25           0 : }
      26             : 
      27          30 : SCMMixingKimAndChung::SCMMixingKimAndChung(const InputParameters & parameters)
      28             :   : SCMMixingClosureBase(parameters),
      29          30 :     _is_tri_lattice(dynamic_cast<const TriSubChannelMesh *>(&_subchannel_mesh) != nullptr),
      30          30 :     _sch_mesh(_subchannel_mesh),
      31          30 :     _S_soln(_subproblem.getVariable(0, "S")),
      32          30 :     _mdot_soln(_subproblem.getVariable(0, "mdot")),
      33          30 :     _w_perim_soln(_subproblem.getVariable(0, "w_perim")),
      34          30 :     _mu_soln(_subproblem.getVariable(0, "mu")),
      35          30 :     _P_soln(_subproblem.getVariable(0, "P")),
      36          60 :     _T_soln(_subproblem.getVariable(0, "T"))
      37             : {
      38          30 :   if (_is_tri_lattice)
      39             :   {
      40          24 :     const auto & tri_mesh = static_cast<const TriSubChannelMesh &>(_subchannel_mesh);
      41          24 :     if (tri_mesh.getWireDiameter() != 0.0 || tri_mesh.getWireLeadLength() != 0.0)
      42           0 :       mooseError("SCMMixingKimAndChung applies only to bare-pin assemblies and cannot be used "
      43             :                  "with wire-wrapped triangular bundles.");
      44             :   }
      45          30 : }
      46             : 
      47             : Real
      48     3390000 : SCMMixingKimAndChung::computeMixingParameter(const unsigned int i_gap, const unsigned int iz) const
      49             : {
      50     3390000 :   if (_is_tri_lattice)
      51      702000 :     return computeTriLatticeMixingParameter(i_gap, iz);
      52             :   else
      53     2688000 :     return computeQuadLatticeMixingParameter(i_gap, iz);
      54             : }
      55             : 
      56             : /**
      57             :  * Common implementation for both tri and quad lattices.
      58             :  * The only lattice-dependent quantities are delta and sf (shape factor).
      59             :  */
      60             : Real
      61     3390000 : SCMMixingKimAndChung::computeLatticeMixingParameter(const unsigned int i_gap,
      62             :                                                     const unsigned int iz,
      63             :                                                     const Real delta,
      64             :                                                     const Real sf) const
      65             : {
      66     3390000 :   const Real P_out = _scm_problem.getOutletPressure();
      67             :   const auto fp = _scm_problem.getSinglePhaseFluidProperties();
      68     3390000 :   const Real pin_diameter = _sch_mesh.getPinDiameter();
      69             : 
      70     3390000 :   const auto chans = _sch_mesh.getGapChannels(i_gap);
      71     3390000 :   const unsigned int i_ch = chans.first;
      72     3390000 :   const unsigned int j_ch = chans.second;
      73             : 
      74     3390000 :   const Node * const node_in_i = _sch_mesh.getChannelNode(i_ch, iz - 1);
      75     3390000 :   const Node * const node_out_i = _sch_mesh.getChannelNode(i_ch, iz);
      76     3390000 :   const Node * const node_in_j = _sch_mesh.getChannelNode(j_ch, iz - 1);
      77     3390000 :   const Node * const node_out_j = _sch_mesh.getChannelNode(j_ch, iz);
      78             : 
      79     3390000 :   const Real Si_in = _S_soln(node_in_i);
      80     3390000 :   const Real Sj_in = _S_soln(node_in_j);
      81     3390000 :   const Real Si_out = _S_soln(node_out_i);
      82     3390000 :   const Real Sj_out = _S_soln(node_out_j);
      83             : 
      84     3390000 :   const Real Si = 0.5 * (Si_in + Si_out);
      85     3390000 :   const Real Sj = 0.5 * (Sj_in + Sj_out);
      86     3390000 :   const Real S_total = Si_in + Sj_in + Si_out + Sj_out;
      87             : 
      88             :   // crossflow area between channels i,j (dz*gap_width)
      89     3390000 :   const Real gap = _sch_mesh.getGapWidth(iz, i_gap);
      90             : 
      91     3390000 :   const Real massflux_i = 0.5 * (_mdot_soln(node_in_i) + _mdot_soln(node_out_i)) / Si;
      92     3390000 :   const Real massflux_j = 0.5 * (_mdot_soln(node_in_j) + _mdot_soln(node_out_j)) / Sj;
      93             : 
      94             :   // average massflux definition
      95             :   const Real avg_massflux =
      96     3390000 :       0.5 * ((_mdot_soln(node_in_i) + _mdot_soln(node_in_j)) / (Si_in + Sj_in) +
      97     3390000 :              (_mdot_soln(node_out_i) + _mdot_soln(node_out_j)) / (Si_out + Sj_out));
      98             : 
      99     3390000 :   const Real w_perim_i = 0.5 * (_w_perim_soln(node_in_i) + _w_perim_soln(node_out_i));
     100     3390000 :   const Real w_perim_j = 0.5 * (_w_perim_soln(node_in_j) + _w_perim_soln(node_out_j));
     101             : 
     102     3390000 :   const Real mu_i = 0.5 * (_mu_soln(node_in_i) + _mu_soln(node_out_i));
     103     3390000 :   const Real mu_j = 0.5 * (_mu_soln(node_in_j) + _mu_soln(node_out_j));
     104     3390000 :   const Real avg_mu = 0.5 * (mu_i + mu_j);
     105             : 
     106     3390000 :   const Real avg_hD = 4.0 * (Si + Sj) / (w_perim_i + w_perim_j);
     107     3390000 :   const Real hD_i = 4.0 * Si / w_perim_i;
     108     3390000 :   const Real hD_j = 4.0 * Sj / w_perim_j;
     109             : 
     110     3390000 :   const Real Re_i = massflux_i * hD_i / mu_i;
     111     3390000 :   const Real Re_j = massflux_j * hD_j / mu_j;
     112     3390000 :   const Real Re = avg_massflux * avg_hD / avg_mu;
     113             : 
     114             :   constexpr Real gamma = 20.0; // empirical constant
     115             :   constexpr Real a = 0.18;
     116             :   constexpr Real b = 0.2;
     117             : 
     118             :   const auto friction_args_i = FrictionStruct(i_ch, Re_i, Si, w_perim_i);
     119             :   const auto friction_args_j = FrictionStruct(j_ch, Re_j, Sj, w_perim_j);
     120             : 
     121     3390000 :   const Real fi = _scm_problem.getFrictionClosure()->computeFrictionFactor(friction_args_i);
     122     3390000 :   const Real fj = _scm_problem.getFrictionClosure()->computeFrictionFactor(friction_args_j);
     123     3390000 :   const Real f = 0.5 * (fi + fj);
     124             : 
     125             :   // average heat conduction
     126             :   const Real avg_k =
     127     3390000 :       (1.0 / S_total) * (fp->k_from_p_T(_P_soln(node_out_i) + P_out, _T_soln(node_out_i)) * Si_out +
     128     3390000 :                          fp->k_from_p_T(_P_soln(node_in_i) + P_out, _T_soln(node_in_i)) * Si_in +
     129     3390000 :                          fp->k_from_p_T(_P_soln(node_out_j) + P_out, _T_soln(node_out_j)) * Sj_out +
     130     3390000 :                          fp->k_from_p_T(_P_soln(node_in_j) + P_out, _T_soln(node_in_j)) * Sj_in);
     131             : 
     132             :   // average specific heat capacity
     133             :   const Real avg_cp = (1.0 / S_total) *
     134     3390000 :                       (fp->cp_from_p_T(_P_soln(node_out_i) + P_out, _T_soln(node_out_i)) * Si_out +
     135     3390000 :                        fp->cp_from_p_T(_P_soln(node_in_i) + P_out, _T_soln(node_in_i)) * Si_in +
     136     3390000 :                        fp->cp_from_p_T(_P_soln(node_out_j) + P_out, _T_soln(node_out_j)) * Sj_out +
     137     3390000 :                        fp->cp_from_p_T(_P_soln(node_in_j) + P_out, _T_soln(node_in_j)) * Sj_in);
     138             : 
     139     3390000 :   const Real Pr = avg_mu * avg_cp / avg_k;                  // Prandtl number
     140     3390000 :   const Real Pr_t = Pr * (Re / gamma) * std::sqrt(f / 8.0); // Turbulent Prandtl number
     141             : 
     142             :   // This form is intended for liquid-metal applications; warn when the local average Prandtl
     143             :   // number is outside the usual low-Pr liquid-metal regime.
     144     3390000 :   if (Pr >= 0.1)
     145     2688004 :     flagSolutionWarning("Prandtl number (Pr) outside the low-Prandtl-number liquid-metal range "
     146             :                         "expected by the Kim-Chung turbulent mixing closure.");
     147             : 
     148     3390000 :   const Real L_x = sf * delta;  // axial length scale (gap is the lateral length scale)
     149     3390000 :   const Real lamda = gap / L_x; // aspect ratio
     150     3390000 :   const Real a_x = 1.0 - 2.0 * lamda * lamda / libMesh::pi; // velocity coefficient
     151             : 
     152             :   const Real z_FP_over_D =
     153     3390000 :       (2.0 * L_x / pin_diameter) *
     154     3390000 :       (1.0 + (-0.5 * std::log(lamda) + 0.5 * std::log(4.0) - 0.25) * lamda * lamda);
     155             : 
     156     3390000 :   const Real Str =
     157     3390000 :       1.0 / (0.822 * (gap / pin_diameter) + 0.144); // Strouhal number (Wu & Trupp 1994)
     158             : 
     159     3390000 :   const Real freq_factor = 2.0 / Utility::pow<2>(gamma) * std::sqrt(a / 8.0) * (avg_hD / gap);
     160             : 
     161     3390000 :   const Real rod_mixing = (1.0 / Pr_t) * lamda;
     162     3390000 :   const Real axial_mixing = a_x * z_FP_over_D * Str;
     163             : 
     164             :   // Mixing Stanton number following Kim and Chung (2001), Eq. 25, as used by
     165             :   // Jeong et al. (2005), Eq. 19:
     166             :   //
     167             :   // St_g = 2/gamma^2 * sqrt(a/8) * (D_h/g)
     168             :   //       * [lambda/Pr_t + a_x * (z_FP/D) * Str] * Re^(-b/2).
     169             :   //
     170             :   // The subchannel problem converts St_g to the turbulent interchange flow through
     171             :   // w'_ij = beta * S_ij * G_bar (St_g = beta). This correlation doesn't include the
     172             :   // molecular diffusion term which represents the conductive heat transfer,
     173             :   // which is separately modeled in the subchannel energy equation.
     174     3390000 :   return freq_factor * (rod_mixing + axial_mixing) * std::pow(Re, -b / 2.0);
     175             : }
     176             : 
     177             : Real
     178      702000 : SCMMixingKimAndChung::computeTriLatticeMixingParameter(const unsigned int i_gap,
     179             :                                                        const unsigned int iz) const
     180             : {
     181             :   // main differences for tri lattice: delta and sf
     182      702000 :   const Real pitch = _sch_mesh.getPitch();
     183      702000 :   const Real delta = pitch / std::sqrt(3.0); // centroid to centroid distance
     184             :   constexpr Real sf = 2.0 / 3.0;             // shape factor
     185      702000 :   return computeLatticeMixingParameter(i_gap, iz, delta, sf);
     186             : }
     187             : 
     188             : Real
     189     2688000 : SCMMixingKimAndChung::computeQuadLatticeMixingParameter(const unsigned int i_gap,
     190             :                                                         const unsigned int iz) const
     191             : {
     192             :   // main differences for quad lattice: delta and sf
     193     2688000 :   const Real pitch = _sch_mesh.getPitch();
     194             :   const Real delta = pitch; // centroid to centroid distance
     195             :   constexpr Real sf = 1.0;  // shape factor
     196     2688000 :   return computeLatticeMixingParameter(i_gap, iz, delta, sf);
     197             : }