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 "SCMFrictionUpdatedChengTodreas.h"
      11             : 
      12             : registerMooseObject("SubChannelApp", SCMFrictionUpdatedChengTodreas);
      13             : 
      14             : InputParameters
      15         359 : SCMFrictionUpdatedChengTodreas::validParams()
      16             : {
      17         359 :   InputParameters params = SCMFrictionClosureBase::validParams();
      18         359 :   params.addClassDescription("Class that computes the axial friction factor using the updated "
      19             :                              "Cheng Todreas correlations.");
      20         359 :   return params;
      21           0 : }
      22             : 
      23         191 : SCMFrictionUpdatedChengTodreas::SCMFrictionUpdatedChengTodreas(const InputParameters & parameters)
      24             :   : SCMFrictionClosureBase(parameters),
      25         191 :     _is_tri_lattice(dynamic_cast<const TriSubChannelMesh *>(&_subchannel_mesh) != nullptr),
      26         191 :     _tri_sch_mesh(dynamic_cast<const TriSubChannelMesh *>(&_subchannel_mesh)),
      27         191 :     _quad_sch_mesh(dynamic_cast<const QuadSubChannelMesh *>(&_subchannel_mesh)),
      28         191 :     _has_wire_wrap(_is_tri_lattice && _tri_sch_mesh->getWireDiameter() != 0.0 &&
      29         349 :                    _tri_sch_mesh->getWireLeadLength() != 0.0)
      30             : {
      31         191 :   if (_is_tri_lattice &&
      32         182 :       ((_tri_sch_mesh->getWireDiameter() != 0.0) != (_tri_sch_mesh->getWireLeadLength() != 0.0)))
      33           0 :     mooseError(name(),
      34             :                ": wire-wrapped bundle friction requires both wire diameter and wire lead length. "
      35             :                "Set both to zero for a bare pin bundle.");
      36             : 
      37         191 :   if (_is_tri_lattice)
      38             :   {
      39         182 :     const auto pitch = _subchannel_mesh.getPitch();
      40         182 :     const auto pin_diameter = _subchannel_mesh.getPinDiameter();
      41         182 :     const auto p_over_d = pitch / pin_diameter;
      42         182 :     const auto wire_lead_to_diameter = _tri_sch_mesh->getWireLeadLength() / pin_diameter;
      43         182 :     const unsigned int Nr = _tri_sch_mesh->getNumOfRings();
      44         182 :     const unsigned int num_pins = 1 + 3 * Nr * (Nr - 1);
      45             : 
      46             :     // The updated Cheng-Todreas detailed triangular friction factor correlation is based on
      47             :     // data spanning 1.0 <= P/D <= 1.42, 4 <= H/D <= 52, 7 <= Npin <= 271,
      48             :     // and 50 <= Re <= 1e6.
      49         182 :     if (p_over_d < 1.0 || p_over_d > 1.42)
      50           0 :       flagSolutionWarning("Pitch-over-pin diameter ratio (P/D) outside the updated "
      51             :                           "Cheng-Todreas friction correlation data range.");
      52         182 :     if (_has_wire_wrap && (wire_lead_to_diameter < 4.0 || wire_lead_to_diameter > 52.0))
      53         132 :       flagSolutionWarning("Wire lead length-over-pin diameter ratio (H/D) outside the updated "
      54             :                           "Cheng-Todreas friction correlation data range.");
      55         182 :     if (num_pins < 7 || num_pins > 271)
      56           0 :       flagSolutionWarning("Number of pins outside the updated Cheng-Todreas friction correlation "
      57             :                           "data range.");
      58             :   }
      59         191 : }
      60             : 
      61             : Real
      62    14201085 : SCMFrictionUpdatedChengTodreas::computeFrictionFactor(const FrictionStruct & friction_args) const
      63             : {
      64    14201085 :   if (_is_tri_lattice)
      65     7175685 :     return computeTriLatticeFrictionFactor(friction_args);
      66             :   else
      67     7025400 :     return computeQuadLatticeFrictionFactor(friction_args);
      68             : }
      69             : 
      70             : Real
      71     7175685 : SCMFrictionUpdatedChengTodreas::computeTriLatticeFrictionFactor(
      72             :     const FrictionStruct & friction_args) const
      73             : {
      74     7175685 :   const auto Re = friction_args.Re;
      75     7175685 :   const auto i_ch = friction_args.i_ch;
      76     7175685 :   const auto S = friction_args.S;
      77     7175685 :   const auto w_perim = friction_args.w_perim;
      78     7175685 :   const auto Dh_i = 4.0 * S / w_perim;
      79             :   Real aL, b1L, b2L, cL;
      80             :   Real aT, b1T, b2T, cT;
      81     7175685 :   const Real & pitch = _subchannel_mesh.getPitch();
      82     7175685 :   const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
      83     7175685 :   const Real & wire_lead_length = _tri_sch_mesh->getWireLeadLength();
      84             :   const Real & wire_diameter = _tri_sch_mesh->getWireDiameter();
      85     7175685 :   const auto p_over_d = pitch / pin_diameter;
      86     7175685 :   const auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
      87             :   // This gap is a constant value for the whole assembly. Might want to make it
      88             :   // subchannel specific in the future if we have duct deformation.
      89     7175685 :   const auto gap = _tri_sch_mesh->getDuctToPinGap();
      90     7175685 :   const auto w_over_d = (pin_diameter + gap) / pin_diameter;
      91             : 
      92     7175685 :   if (Re < 50.0 || Re > 1.0e6)
      93       14367 :     flagSolutionWarning("Reynolds number (Re) outside the updated Cheng-Todreas friction "
      94             :                         "correlation data range.");
      95             : 
      96     7175685 :   const auto ReL = std::pow(10, (p_over_d - 1)) * 320.0;
      97     7175685 :   const auto ReT = std::pow(10, 0.7 * (p_over_d - 1)) * 1.0E+4;
      98     7175685 :   const auto psi = std::log(Re / ReL) / std::log(ReT / ReL);
      99             : 
     100             :   // Find the coefficients of bare Pin bundle friction factor
     101             :   // correlations for turbulent and laminar flow regimes. Todreas & Kazimi, Nuclear Systems
     102             :   // second edition, Volume 1, Chapter 9.6
     103     7175685 :   if (subch_type == EChannelType::CENTER)
     104             :   {
     105     4608591 :     if (p_over_d < 1.1)
     106             :     {
     107             :       aL = 26.0;
     108             :       b1L = 888.2;
     109             :       b2L = -3334.0;
     110             :       aT = 0.09378;
     111             :       b1T = 1.398;
     112             :       b2T = -8.664;
     113             :     }
     114             :     else
     115             :     {
     116             :       aL = 62.97;
     117             :       b1L = 216.9;
     118             :       b2L = -190.2;
     119             :       aT = 0.1458;
     120             :       b1T = 0.03632;
     121             :       b2T = -0.03333;
     122             :     }
     123             :     // laminar flow friction factor for bare Pin bundle - Center subchannel
     124     4608591 :     cL = aL + b1L * (p_over_d - 1) + b2L * Utility::pow<2>((p_over_d - 1));
     125             :     // turbulent flow friction factor for bare Pin bundle - Center subchannel
     126     4608591 :     cT = aT + b1T * (p_over_d - 1) + b2T * Utility::pow<2>((p_over_d - 1));
     127             :   }
     128     2567094 :   else if (subch_type == EChannelType::EDGE)
     129             :   {
     130     1827498 :     if (w_over_d < 1.1)
     131             :     {
     132             :       aL = 26.18;
     133             :       b1L = 554.5;
     134             :       b2L = -1480.0;
     135             :       aT = 0.09377;
     136             :       b1T = 0.8732;
     137             :       b2T = -3.341;
     138             :     }
     139             :     else
     140             :     {
     141             :       aL = 44.4;
     142             :       b1L = 256.7;
     143             :       b2L = -267.6;
     144             :       aT = 0.1430;
     145             :       b1T = 0.04199;
     146             :       b2T = -0.04428;
     147             :     }
     148             :     // laminar flow friction factor for bare Pin bundle - Edge subchannel
     149     1827498 :     cL = aL + b1L * (w_over_d - 1) + b2L * Utility::pow<2>((w_over_d - 1));
     150             :     // turbulent flow friction factor for bare Pin bundle - Edge subchannel
     151     1827498 :     cT = aT + b1T * (w_over_d - 1) + b2T * Utility::pow<2>((w_over_d - 1));
     152             :   }
     153             :   else
     154             :   {
     155      739596 :     if (w_over_d < 1.1)
     156             :     {
     157             :       aL = 26.98;
     158             :       b1L = 1636.0;
     159             :       b2L = -10050.0;
     160             :       aT = 0.1004;
     161             :       b1T = 1.625;
     162             :       b2T = -11.85;
     163             :     }
     164             :     else
     165             :     {
     166             :       aL = 87.26;
     167             :       b1L = 38.59;
     168             :       b2L = -55.12;
     169             :       aT = 0.1499;
     170             :       b1T = 0.006706;
     171             :       b2T = -0.009567;
     172             :     }
     173             :     // laminar flow friction factor for bare Pin bundle - Corner subchannel
     174      739596 :     cL = aL + b1L * (w_over_d - 1) + b2L * Utility::pow<2>((w_over_d - 1));
     175             :     // turbulent flow friction factor for bare Pin bundle - Corner subchannel
     176      739596 :     cT = aT + b1T * (w_over_d - 1) + b2T * Utility::pow<2>((w_over_d - 1));
     177             :   }
     178             : 
     179             :   // Find the coefficients of wire-wrapped Pin bundle friction factor
     180             :   // correlations for turbulent and laminar flow regimes. Todreas & Kazimi, Nuclear Systems
     181             :   // Volume 1 Chapter 9-6 also Chen and Todreas (2018).
     182     7175685 :   if (_has_wire_wrap)
     183             :   {
     184             :     const auto theta =
     185     5417925 :         std::acos(wire_lead_length /
     186     5417925 :                   std::sqrt(Utility::pow<2>(wire_lead_length) +
     187     5417925 :                             Utility::pow<2>(libMesh::pi * (pin_diameter + wire_diameter))));
     188     5417925 :     const auto wd_t = (19.56 - 98.71 * (wire_diameter / pin_diameter) +
     189     5417925 :                        303.47 * Utility::pow<2>((wire_diameter / pin_diameter))) *
     190     5417925 :                       std::pow((wire_lead_length / pin_diameter), -0.541);
     191     5417925 :     const auto wd_l = 1.4 * wd_t;
     192     5417925 :     const auto ws_t = -11.0 * std::log(wire_lead_length / pin_diameter) + 19.0;
     193             :     const auto ws_l = ws_t;
     194             :     Real ar = 0.0;
     195             :     Real a_p = 0.0;
     196             : 
     197     5417925 :     if (subch_type == EChannelType::CENTER)
     198             :     {
     199             :       // wetted perimeter for center subchannel and bare Pin bundle
     200     3465771 :       const Real pw_p = libMesh::pi * pin_diameter / 2.0;
     201             :       // wire projected area - center subchannel wire-wrapped bundle
     202     3465771 :       ar = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
     203             :       // bare Pin bundle center subchannel flow area (normal area + wire area)
     204     3465771 :       a_p = S + libMesh::pi * Utility::pow<2>(wire_diameter) / 8.0 / std::cos(theta);
     205             :       // turbulent friction factor equation constant - Center subchannel
     206     3465771 :       cT *= (pw_p / w_perim);
     207     3465771 :       cT += wd_t * (3.0 * ar / a_p) * (Dh_i / wire_lead_length) *
     208     3465771 :             std::pow((Dh_i / wire_diameter), 0.18);
     209             :       // laminar friction factor equation constant - Center subchannel
     210     3465771 :       cL *= (pw_p / w_perim);
     211     3465771 :       cL += wd_l * (3.0 * ar / a_p) * (Dh_i / wire_lead_length) * (Dh_i / wire_diameter);
     212             :     }
     213     1952154 :     else if (subch_type == EChannelType::EDGE)
     214             :     {
     215             :       // wire projected area - edge subchannel wire-wrapped bundle
     216     1366548 :       ar = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 4.0;
     217             :       // bare Pin bundle edge subchannel flow area (normal area + wire area)
     218     1366548 :       a_p = S + libMesh::pi * Utility::pow<2>(wire_diameter) / 8.0 / std::cos(theta);
     219             :       // turbulent friction factor equation constant - Edge subchannel
     220     1366548 :       cT *= std::pow((1 + ws_t * (ar / a_p) * Utility::pow<2>(std::tan(theta))), 1.41);
     221             :       // laminar friction factor equation constant - Edge subchannel
     222     1366548 :       cL *= (1 + ws_l * (ar / a_p) * Utility::pow<2>(std::tan(theta)));
     223             :     }
     224             :     else
     225             :     {
     226             :       // wire projected area - corner subchannel wire-wrapped bundle
     227      585606 :       ar = libMesh::pi * (pin_diameter + wire_diameter) * wire_diameter / 6.0;
     228             :       // bare Pin bundle corner subchannel flow area (normal area + wire area)
     229      585606 :       a_p = S + libMesh::pi * Utility::pow<2>(wire_diameter) / 24.0 / std::cos(theta);
     230             :       // turbulent friction factor equation constant - Corner subchannel
     231      585606 :       cT *= std::pow((1 + ws_t * (ar / a_p) * Utility::pow<2>(std::tan(theta))), 1.41);
     232             :       // laminar friction factor equation constant - Corner subchannel
     233      585606 :       cL *= (1 + ws_l * (ar / a_p) * Utility::pow<2>(std::tan(theta)));
     234             :     }
     235             :   }
     236             :   // laminar friction factor and turbulent friction factor coefficients
     237             :   const Real bL = -1.0;
     238             :   const Real bT = -0.18;
     239     7175685 :   auto fL = cL * std::pow(Re, bL);
     240     7175685 :   auto fT = cT * std::pow(Re, bT);
     241             : 
     242     7175685 :   if (Re < ReL)
     243             :   {
     244             :     // laminar flow
     245             :     return fL;
     246             :   }
     247     7161321 :   else if (Re > ReT)
     248             :   {
     249             :     // turbulent flow
     250             :     return fT;
     251             :   }
     252             :   else
     253             :   {
     254             :     // transient flow: psi definition uses a Bulk ReT/ReL number, same for all channels
     255     1127535 :     return fL * std::pow((1 - psi), 1.0 / 3.0) * (1 - std::pow(psi, 7)) +
     256     1127535 :            fT * std::pow(psi, 1.0 / 3.0);
     257             :   }
     258             : }
     259             : 
     260             : Real
     261     7025400 : SCMFrictionUpdatedChengTodreas::computeQuadLatticeFrictionFactor(
     262             :     const FrictionStruct & friction_args) const
     263             : {
     264     7025400 :   const auto Re = friction_args.Re;
     265     7025400 :   const auto i_ch = friction_args.i_ch;
     266             :   /// Todreas-Kazimi NUCLEAR SYSTEMS, second edition, Volume 1, 2011
     267             :   Real aL, b1L, b2L, cL;
     268             :   Real aT, b1T, b2T, cT;
     269     7025400 :   const auto pitch = _subchannel_mesh.getPitch();
     270     7025400 :   const auto pin_diameter = _subchannel_mesh.getPinDiameter();
     271             :   // This gap is a constant value for the whole assembly. Might want to make it
     272             :   // subchannel specific in the future if we have duct deformation.
     273     7025400 :   const auto side_gap = _quad_sch_mesh->getSideGap();
     274     7025400 :   const auto w = (pin_diameter / 2.0) + (pitch / 2.0) + side_gap;
     275     7025400 :   const auto p_over_d = pitch / pin_diameter;
     276     7025400 :   const auto w_over_d = w / pin_diameter;
     277     7025400 :   const auto ReL = std::pow(10, (p_over_d - 1)) * 320.0;
     278     7025400 :   const auto ReT = std::pow(10, 0.7 * (p_over_d - 1)) * 1.0E+4;
     279     7025400 :   const auto psi = std::log(Re / ReL) / std::log(ReT / ReL);
     280     7025400 :   const auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
     281             : 
     282             :   // Find the coefficients of bare Pin bundle friction factor
     283             :   // correlations for turbulent and laminar flow regimes. Todreas & Kazimi, Nuclear Systems Volume
     284             :   // 1
     285     7025400 :   if (subch_type == EChannelType::CENTER)
     286             :   {
     287     3620600 :     if (p_over_d < 1.1)
     288             :     {
     289             :       aL = 26.37;
     290             :       b1L = 374.2;
     291             :       b2L = -493.9;
     292             :       aT = 0.09423;
     293             :       b1T = 0.5806;
     294             :       b2T = -1.239;
     295             :     }
     296             :     else
     297             :     {
     298             :       aL = 35.55;
     299             :       b1L = 263.7;
     300             :       b2L = -190.2;
     301             :       aT = 0.1339;
     302             :       b1T = 0.09059;
     303             :       b2T = -0.09926;
     304             :     }
     305             :     // laminar flow friction factor for bare Pin bundle - Center subchannel
     306     3620600 :     cL = aL + b1L * (p_over_d - 1) + b2L * Utility::pow<2>((p_over_d - 1));
     307             :     // turbulent flow friction factor for bare Pin bundle - Center subchannel
     308     3620600 :     cT = aT + b1T * (p_over_d - 1) + b2T * Utility::pow<2>((p_over_d - 1));
     309             :   }
     310     3404800 :   else if (subch_type == EChannelType::EDGE)
     311             :   {
     312     2867200 :     if (p_over_d < 1.1)
     313             :     {
     314             :       aL = 26.18;
     315             :       b1L = 554.5;
     316             :       b2L = -1480;
     317             :       aT = 0.09377;
     318             :       b1T = 0.8732;
     319             :       b2T = -3.341;
     320             :     }
     321             :     else
     322             :     {
     323             :       aL = 44.40;
     324             :       b1L = 256.7;
     325             :       b2L = -267.6;
     326             :       aT = 0.1430;
     327             :       b1T = 0.04199;
     328             :       b2T = -0.04428;
     329             :     }
     330             :     // laminar flow friction factor for bare Pin bundle - Edge subchannel
     331     2867200 :     cL = aL + b1L * (w_over_d - 1) + b2L * Utility::pow<2>((w_over_d - 1));
     332             :     // turbulent flow friction factor for bare Pin bundle - Edge subchannel
     333     2867200 :     cT = aT + b1T * (w_over_d - 1) + b2T * Utility::pow<2>((w_over_d - 1));
     334             :   }
     335             :   else
     336             :   {
     337      537600 :     if (p_over_d < 1.1)
     338             :     {
     339             :       aL = 28.62;
     340             :       b1L = 715.9;
     341             :       b2L = -2807;
     342             :       aT = 0.09755;
     343             :       b1T = 1.127;
     344             :       b2T = -6.304;
     345             :     }
     346             :     else
     347             :     {
     348             :       aL = 58.83;
     349             :       b1L = 160.7;
     350             :       b2L = -203.5;
     351             :       aT = 0.1452;
     352             :       b1T = 0.02681;
     353             :       b2T = -0.03411;
     354             :     }
     355             :     // laminar flow friction factor for bare Pin bundle - Corner subchannel
     356      537600 :     cL = aL + b1L * (w_over_d - 1) + b2L * Utility::pow<2>((w_over_d - 1));
     357             :     // turbulent flow friction factor for bare Pin bundle - Corner subchannel
     358      537600 :     cT = aT + b1T * (w_over_d - 1) + b2T * Utility::pow<2>((w_over_d - 1));
     359             :   }
     360             :   // laminar friction factor and turbulent friction factor coefficients
     361             :   const Real bL = -1.0;
     362             :   const Real bT = -0.18;
     363     7025400 :   auto fL = cL * std::pow(Re, bL);
     364     7025400 :   auto fT = cT * std::pow(Re, bT);
     365             : 
     366     7025400 :   if (Re < ReL)
     367             :   {
     368             :     // laminar flow
     369             :     return fL;
     370             :   }
     371     7025400 :   else if (Re > ReT)
     372             :   {
     373             :     // turbulent flow
     374             :     return fT;
     375             :   }
     376             :   else
     377             :   {
     378             :     // transient flow: psi definition uses a Bulk ReT/ReL number, same for all channels
     379           0 :     return fL * std::pow((1 - psi), 1.0 / 3.0) * (1 - std::pow(psi, 7)) +
     380           0 :            fT * std::pow(psi, 1.0 / 3.0);
     381             :   }
     382             : }