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 "SCMMixingChengTodreas.h"
      11             : 
      12             : registerMooseObject("SubChannelApp", SCMMixingChengTodreas);
      13             : 
      14             : InputParameters
      15         258 : SCMMixingChengTodreas::validParams()
      16             : {
      17         258 :   InputParameters params = SCMMixingClosureBase::validParams();
      18         258 :   params.addClassDescription("Class that models the turbulent mixing coefficient for wire-wrapped "
      19             :                              "triangular assemblies using the Cheng Todreas correlations.");
      20         258 :   return params;
      21           0 : }
      22             : 
      23         138 : SCMMixingChengTodreas::SCMMixingChengTodreas(const InputParameters & parameters)
      24             :   : SCMMixingClosureBase(parameters),
      25         138 :     _is_tri_lattice(dynamic_cast<const TriSubChannelMesh *>(&_subchannel_mesh) != nullptr),
      26         138 :     _tri_sch_mesh(dynamic_cast<const TriSubChannelMesh *>(&_subchannel_mesh)),
      27         138 :     _S_soln(_subproblem.getVariable(0, "S")),
      28         138 :     _mdot_soln(_subproblem.getVariable(0, "mdot")),
      29         138 :     _w_perim_soln(_subproblem.getVariable(0, "w_perim")),
      30         276 :     _mu_soln(_subproblem.getVariable(0, "mu"))
      31             : {
      32         138 :   if (!_is_tri_lattice)
      33           0 :     mooseError("This correlation applies only for triangular assemblies");
      34             : 
      35         138 :   if (_tri_sch_mesh->getWireLeadLength() == 0 || _tri_sch_mesh->getWireDiameter() == 0)
      36           0 :     mooseError("This correlation applies only for wire-wrapped assemblies");
      37             : 
      38         138 :   const Real pitch = _subchannel_mesh.getPitch();
      39         138 :   const Real pin_diameter = _subchannel_mesh.getPinDiameter();
      40         138 :   const Real pitch_to_diameter = pitch / pin_diameter;
      41         138 :   const Real wire_lead_to_diameter = _tri_sch_mesh->getWireLeadLength() / pin_diameter;
      42         138 :   const unsigned int Nr = _tri_sch_mesh->getNumOfRings();
      43         138 :   const unsigned int num_pins = 1 + 3 * Nr * (Nr - 1);
      44             : 
      45             :   // Cheng and Todreas (1986) fitted the wire-wrapped mixing parameters over the ranges
      46             :   // 1.067 <= P/D <= 1.35, 4 <= H/D <= 52, and 7 <= Npin <= 217.
      47         138 :   if (pitch_to_diameter < 1.067 || pitch_to_diameter > 1.35)
      48           0 :     flagSolutionWarning("Pitch-over-pin diameter ratio (P/D) outside the Cheng-Todreas "
      49             :                         "wire-wrapped mixing correlation data range.");
      50         138 :   if (wire_lead_to_diameter < 4.0 || wire_lead_to_diameter > 52.0)
      51         102 :     flagSolutionWarning("Wire lead length-over-pin diameter ratio (H/D) outside the "
      52             :                         "Cheng-Todreas wire-wrapped mixing correlation data range.");
      53         138 :   if (num_pins < 7 || num_pins > 217)
      54           0 :     flagSolutionWarning("Number of pins outside the Cheng-Todreas wire-wrapped mixing "
      55             :                         "correlation data range.");
      56         138 : }
      57             : 
      58             : Real
      59    10052160 : SCMMixingChengTodreas::computeMixingParameter(const unsigned int i_gap, const unsigned int iz) const
      60             : {
      61             :   Real beta = 0.0;
      62             : 
      63    10052160 :   const Real pitch = _subchannel_mesh.getPitch();
      64    10052160 :   const Real pin_diameter = _subchannel_mesh.getPinDiameter();
      65             : 
      66    10052160 :   const Real wire_lead_length = _tri_sch_mesh->getWireLeadLength();
      67    10052160 :   const Real wire_diameter = _tri_sch_mesh->getWireDiameter();
      68    10052160 :   const unsigned int Nr = _tri_sch_mesh->getNumOfRings();
      69             : 
      70    10052160 :   const auto chans = _subchannel_mesh.getGapChannels(i_gap);
      71    10052160 :   const unsigned int i_ch = chans.first;
      72    10052160 :   const unsigned int j_ch = chans.second;
      73             : 
      74    10052160 :   const auto subch_type_i = _subchannel_mesh.getSubchannelType(i_ch);
      75    10052160 :   const auto subch_type_j = _subchannel_mesh.getSubchannelType(j_ch);
      76             : 
      77    10052160 :   const Node * const node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
      78    10052160 :   const Node * const node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
      79    10052160 :   const Node * const node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
      80    10052160 :   const Node * const node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
      81             : 
      82    10052160 :   const Real Si_in = _S_soln(node_in_i);
      83    10052160 :   const Real Sj_in = _S_soln(node_in_j);
      84    10052160 :   const Real Si_out = _S_soln(node_out_i);
      85    10052160 :   const Real Sj_out = _S_soln(node_out_j);
      86             : 
      87    10052160 :   const Real S_total = Si_in + Sj_in + Si_out + Sj_out;
      88    10052160 :   const Real Si = 0.5 * (Si_in + Si_out);
      89    10052160 :   const Real Sj = 0.5 * (Sj_in + Sj_out);
      90             : 
      91    10052160 :   const Real w_perim_i = 0.5 * (_w_perim_soln(node_in_i) + _w_perim_soln(node_out_i));
      92    10052160 :   const Real w_perim_j = 0.5 * (_w_perim_soln(node_in_j) + _w_perim_soln(node_out_j));
      93             : 
      94             :   const Real avg_mu =
      95    10052160 :       (1.0 / S_total) * (_mu_soln(node_out_i) * Si_out + _mu_soln(node_in_i) * Si_in +
      96    10052160 :                          _mu_soln(node_out_j) * Sj_out + _mu_soln(node_in_j) * Sj_in);
      97             : 
      98    10052160 :   const Real avg_hD = 4.0 * (Si + Sj) / (w_perim_i + w_perim_j);
      99             : 
     100             :   const Real avg_massflux =
     101    10052160 :       0.5 * ((_mdot_soln(node_in_i) + _mdot_soln(node_in_j)) / (Si_in + Sj_in) +
     102    10052160 :              (_mdot_soln(node_out_i) + _mdot_soln(node_out_j)) / (Si_out + Sj_out));
     103             : 
     104    10052160 :   const Real Re = avg_massflux * avg_hD / avg_mu;
     105    10052160 :   if (Re < 400.0 || Re > 1.0e6)
     106       53013 :     flagSolutionWarning("Reynolds number (Re) outside the Cheng-Todreas wire-wrapped mixing "
     107             :                         "correlation data range.");
     108             : 
     109             :   // Calculation of flow regime
     110    10052160 :   const Real ReL = 320.0 * std::pow(10.0, pitch / pin_diameter - 1.0);
     111    10052160 :   const Real ReT = 10000.0 * std::pow(10.0, 0.7 * (pitch / pin_diameter - 1.0));
     112             : 
     113             :   // This beta is used by the global turbulent crossflow relation:
     114             :   // w'_ij = beta * S_ij * G_bar. Peripheral sweep flow is handled separately by
     115             :   // computeSweepFlowMixingParameter().
     116    10052160 :   if (subch_type_i == EChannelType::CENTER || subch_type_j == EChannelType::CENTER)
     117             :   {
     118             :     // Calculation of geometric parameters
     119             :     // wire angle
     120             :     const Real theta =
     121     7712280 :         std::acos(wire_lead_length /
     122     7712280 :                   std::sqrt(Utility::pow<2>(wire_lead_length) +
     123     7712280 :                             Utility::pow<2>(libMesh::pi * (pin_diameter + wire_diameter))));
     124             : 
     125             :     // projected area of wire on subchannel
     126     7712280 :     const Real Ar1 = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
     127             : 
     128             :     // bare subchannel flow area
     129     7712280 :     const Real A1prime = (std::sqrt(3.0) / 4.0) * Utility::pow<2>(pitch) -
     130     7712280 :                          libMesh::pi * Utility::pow<2>(pin_diameter) / 8.0;
     131             : 
     132             :     // empirical constant for mixing parameter
     133             :     Real Cm = 0.0;
     134             :     Real CmL_constant = 0.0;
     135             :     Real CmT_constant = 0.0;
     136             : 
     137     7712280 :     if (Nr == 1)
     138             :     {
     139             :       CmT_constant = 0.1;
     140             :       CmL_constant = 0.055;
     141             :     }
     142             :     else
     143             :     {
     144             :       CmT_constant = 0.14;
     145             :       CmL_constant = 0.077;
     146             :     }
     147             : 
     148     7712280 :     const Real CmT = CmT_constant * std::pow((pitch - pin_diameter) / pin_diameter, -0.5);
     149     7712280 :     const Real CmL = CmL_constant * std::pow((pitch - pin_diameter) / pin_diameter, -0.5);
     150             : 
     151     7712280 :     if (Re < ReL)
     152             :     {
     153             :       Cm = CmL;
     154             :     }
     155     7667820 :     else if (Re > ReT)
     156             :     {
     157             :       Cm = CmT;
     158             :     }
     159             :     else
     160             :     {
     161             :       // Simplified intermittency factor; see SCMMixingChengTodreas.md.
     162             :       // Cheng and Todreas (1986) use a more detailed expression for psi.
     163     1537751 :       const Real psi = (std::log(Re) - std::log(ReL)) / (std::log(ReT) - std::log(ReL));
     164             :       const Real gamma = 2.0 / 3.0;
     165     1537751 :       Cm = CmL + (CmT - CmL) * std::pow(psi, gamma);
     166             :     }
     167             : 
     168             :     // mixing parameter
     169     7712280 :     beta = Cm * std::sqrt(Ar1 / A1prime) * std::tan(theta);
     170             :   }
     171    10052160 :   return beta;
     172             : }
     173             : 
     174             : Real
     175     1765680 : SCMMixingChengTodreas::computeSweepFlowMixingParameter(const unsigned int i_gap,
     176             :                                                        const unsigned int iz) const
     177             : {
     178             :   Real beta = 0.0;
     179             : 
     180     1765680 :   const Real pitch = _subchannel_mesh.getPitch();
     181     1765680 :   const Real pin_diameter = _subchannel_mesh.getPinDiameter();
     182             : 
     183     1765680 :   const Real wire_lead_length = _tri_sch_mesh->getWireLeadLength();
     184     1765680 :   const Real wire_diameter = _tri_sch_mesh->getWireDiameter();
     185     1765680 :   const unsigned int Nr = _tri_sch_mesh->getNumOfRings();
     186             : 
     187     1765680 :   const auto chans = _subchannel_mesh.getGapChannels(i_gap);
     188     1765680 :   const unsigned int i_ch = chans.first;
     189     1765680 :   const unsigned int j_ch = chans.second;
     190             : 
     191     1765680 :   const auto subch_type_i = _subchannel_mesh.getSubchannelType(i_ch);
     192     1765680 :   const auto subch_type_j = _subchannel_mesh.getSubchannelType(j_ch);
     193             : 
     194     1765680 :   const Node * const node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
     195     1765680 :   const Node * const node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
     196     1765680 :   const Node * const node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
     197     1765680 :   const Node * const node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
     198             : 
     199     1765680 :   const Real Si_in = _S_soln(node_in_i);
     200     1765680 :   const Real Sj_in = _S_soln(node_in_j);
     201     1765680 :   const Real Si_out = _S_soln(node_out_i);
     202     1765680 :   const Real Sj_out = _S_soln(node_out_j);
     203             : 
     204     1765680 :   const Real S_total = Si_in + Sj_in + Si_out + Sj_out;
     205     1765680 :   const Real Si = 0.5 * (Si_in + Si_out);
     206     1765680 :   const Real Sj = 0.5 * (Sj_in + Sj_out);
     207             : 
     208     1765680 :   const Real w_perim_i = 0.5 * (_w_perim_soln(node_in_i) + _w_perim_soln(node_out_i));
     209     1765680 :   const Real w_perim_j = 0.5 * (_w_perim_soln(node_in_j) + _w_perim_soln(node_out_j));
     210             : 
     211             :   const Real avg_mu =
     212     1765680 :       (1.0 / S_total) * (_mu_soln(node_out_i) * Si_out + _mu_soln(node_in_i) * Si_in +
     213     1765680 :                          _mu_soln(node_out_j) * Sj_out + _mu_soln(node_in_j) * Sj_in);
     214             : 
     215     1765680 :   const Real avg_hD = 4.0 * (Si + Sj) / (w_perim_i + w_perim_j);
     216             : 
     217             :   const Real avg_massflux =
     218     1765680 :       0.5 * ((_mdot_soln(node_in_i) + _mdot_soln(node_in_j)) / (Si_in + Sj_in) +
     219     1765680 :              (_mdot_soln(node_out_i) + _mdot_soln(node_out_j)) / (Si_out + Sj_out));
     220             : 
     221     1765680 :   const Real Re = avg_massflux * avg_hD / avg_mu;
     222     1765680 :   if (Re < 400.0 || Re > 1.0e6)
     223       10263 :     flagSolutionWarning("Reynolds number (Re) outside the Cheng-Todreas wire-wrapped mixing "
     224             :                         "correlation data range.");
     225             : 
     226             :   // Calculation of flow regime
     227     1765680 :   const Real ReL = 320.0 * std::pow(10.0, pitch / pin_diameter - 1.0);
     228     1765680 :   const Real ReT = 10000.0 * std::pow(10.0, 0.7 * (pitch / pin_diameter - 1.0));
     229             : 
     230     1765680 :   if ((subch_type_i == EChannelType::CORNER || subch_type_i == EChannelType::EDGE) &&
     231     1765680 :       (subch_type_j == EChannelType::CORNER || subch_type_j == EChannelType::EDGE))
     232             :   {
     233             :     const Real theta =
     234     1765680 :         std::acos(wire_lead_length /
     235     1765680 :                   std::sqrt(Utility::pow<2>(wire_lead_length) +
     236     1765680 :                             Utility::pow<2>(libMesh::pi * (pin_diameter + wire_diameter))));
     237             : 
     238             :     // Calculation of geometric parameters
     239             :     // distance from pin surface to duct
     240     1765680 :     const Real dpgap = _tri_sch_mesh->getDuctToPinGap();
     241             : 
     242             :     // Edge pitch parameter defined as pin diameter plus distance to duct wall
     243     1765680 :     const Real w = pin_diameter + dpgap;
     244             : 
     245     1765680 :     const Real Ar2 = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 4.0;
     246             : 
     247             :     const Real A2prime =
     248     1765680 :         pitch * (w - pin_diameter / 2.0) - libMesh::pi * Utility::pow<2>(pin_diameter) / 8.0;
     249             : 
     250             :     // empirical constant for mixing parameter
     251             :     Real Cs = 0.0;
     252             :     Real CsL_constant = 0.0;
     253             :     Real CsT_constant = 0.0;
     254             : 
     255     1765680 :     if (Nr == 1)
     256             :     {
     257             :       CsT_constant = 0.6;
     258             :       CsL_constant = 0.33;
     259             :     }
     260             :     else
     261             :     {
     262             :       CsT_constant = 0.75;
     263             :       CsL_constant = 0.413;
     264             :     }
     265             : 
     266     1765680 :     const Real CsL = CsL_constant * std::pow(wire_lead_length / pin_diameter, 0.3);
     267     1765680 :     const Real CsT = CsT_constant * std::pow(wire_lead_length / pin_diameter, 0.3);
     268             : 
     269     1765680 :     if (Re < ReL)
     270             :     {
     271             :       Cs = CsL;
     272             :     }
     273     1755420 :     else if (Re > ReT)
     274             :     {
     275             :       Cs = CsT;
     276             :     }
     277             :     else
     278             :     {
     279             :       // Simplified intermittency factor; see SCMMixingChengTodreas.md.
     280             :       // Cheng and Todreas (1986) use a more detailed expression for psi.
     281      399096 :       const Real psi = (std::log(Re) - std::log(ReL)) / (std::log(ReT) - std::log(ReL));
     282             :       const Real gamma = 2.0 / 3.0;
     283      399096 :       Cs = CsL + (CsT - CsL) * std::pow(psi, gamma);
     284             :     }
     285             : 
     286             :     // Sweep-flow coefficient used only by the peripheral enthalpy calculation.
     287     1765680 :     beta = Cs * std::sqrt(Ar2 / A2prime) * std::tan(theta);
     288             :   }
     289             : 
     290     1765680 :   return beta;
     291             : }