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