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 "SCMTriPowerIC.h"
      11             : #include "Function.h"
      12             : #include "TriSubChannelMesh.h"
      13             : 
      14             : registerMooseObject("SubChannelApp", SCMTriPowerIC);
      15             : registerMooseObjectRenamed("SubChannelApp", TriPowerIC, "06/30/2025 24:00", SCMTriPowerIC);
      16             : 
      17             : InputParameters
      18         239 : SCMTriPowerIC::validParams()
      19             : {
      20         239 :   InputParameters params = TriSubChannelBaseIC::validParams();
      21         239 :   params.addClassDescription(
      22             :       "Computes axial power rate (W/m) that goes into the subchannel cells "
      23             :       "or is assigned to the fuel pins, in a triangular lattice arrangement");
      24         478 :   params.addRequiredParam<PostprocessorName>(
      25             :       "power", "The postprocessor or Real to use for the total power of the subassembly [W]");
      26         478 :   params.addRequiredParam<std::string>(
      27             :       "filename", "name of radial power profile .txt file (should be a single column) [UnitLess].");
      28         478 :   params.addParam<FunctionName>("axial_heat_rate",
      29             :                                 "1.0",
      30             :                                 "user provided normalized function of axial heat rate [Unitless]. "
      31             :                                 "The integral over pin length should equal the heated length");
      32         239 :   return params;
      33           0 : }
      34             : 
      35         134 : SCMTriPowerIC::SCMTriPowerIC(const InputParameters & params)
      36             :   : TriSubChannelBaseIC(params),
      37         134 :     _power(getPostprocessorValue("power")),
      38         134 :     _numberoflines(0),
      39         402 :     _filename(getParam<std::string>("filename")),
      40         268 :     _axial_heat_rate(getFunction("axial_heat_rate"))
      41             : {
      42         134 :   auto n_pins = _mesh.getNumOfPins();
      43         134 :   auto heated_length = _mesh.getHeatedLength();
      44             : 
      45         134 :   _power_dis.resize(n_pins, 1);
      46         134 :   _power_dis.setZero();
      47         134 :   _pin_power_correction.resize(n_pins, 1);
      48         134 :   _pin_power_correction.setOnes();
      49             : 
      50             :   Real vin;
      51         134 :   std::ifstream inFile;
      52             : 
      53         134 :   inFile.open(_filename);
      54         134 :   if (!inFile)
      55           0 :     mooseError(name(), ": Unable to open file: ", _filename);
      56             : 
      57        4744 :   while (inFile >> vin)
      58        4610 :     _numberoflines += 1;
      59             : 
      60         134 :   if (inFile.fail() && !inFile.eof())
      61           0 :     mooseError(name(), ": Non numerical input at line: ", _numberoflines);
      62             : 
      63         134 :   if (_numberoflines != n_pins)
      64           0 :     mooseError(name(), ": Radial profile file doesn't have correct size: ", n_pins);
      65         134 :   inFile.close();
      66             : 
      67         134 :   inFile.open(_filename);
      68             :   int i = 0;
      69        4744 :   while (inFile >> vin)
      70             :   {
      71        4610 :     _power_dis(i, 0) = vin;
      72        4610 :     i++;
      73             :   }
      74         134 :   inFile.close();
      75             : 
      76         134 :   auto sum = _power_dis.sum();
      77             :   // full (100%) power of one pin [W]
      78         134 :   auto fpin_power = _power / sum;
      79             :   // actual pin power [W]
      80         134 :   _ref_power = _power_dis * fpin_power;
      81             :   // Convert the actual pin power to a linear heat rate [W/m]
      82         134 :   _ref_qprime = _ref_power / heated_length;
      83         134 : }
      84             : 
      85             : void
      86         134 : SCMTriPowerIC::initialSetup()
      87             : {
      88         134 :   auto n_pins = _mesh.getNumOfPins();
      89         134 :   auto nz = _mesh.getNumOfAxialCells();
      90         134 :   auto z_grid = _mesh.getZGrid();
      91         134 :   auto heated_length = _mesh.getHeatedLength();
      92         134 :   auto unheated_length_entry = _mesh.getHeatedLengthEntry();
      93             : 
      94         134 :   _estimate_power.resize(n_pins, 1);
      95         134 :   _estimate_power.setZero();
      96        3778 :   for (unsigned int iz = 1; iz < nz + 1; iz++)
      97             :   {
      98             :     // Compute axial location of nodes.
      99        3644 :     auto z2 = z_grid[iz];
     100        3644 :     auto z1 = z_grid[iz - 1];
     101        3644 :     Point p1(0, 0, z1 - unheated_length_entry);
     102        3644 :     Point p2(0, 0, z2 - unheated_length_entry);
     103        3644 :     auto heat1 = _axial_heat_rate.value(_t, p1);
     104        3644 :     auto heat2 = _axial_heat_rate.value(_t, p2);
     105        3644 :     if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
     106        2808 :         MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
     107             :     {
     108             :       // cycle through pins
     109       68148 :       for (unsigned int i_pin = 0; i_pin < n_pins; i_pin++)
     110             :       {
     111             :         // Compute the height of this element.
     112       65940 :         auto dz = z2 - z1;
     113             : 
     114             :         // calculation of power for the first heated segment if nodes don't align
     115       65940 :         if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
     116             :             MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry))
     117             :         {
     118             :           heat1 = 0.0;
     119             :         }
     120             : 
     121             :         // calculation of power for the last heated segment if nodes don't align
     122       65940 :         if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry + heated_length) &&
     123             :             MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
     124             :         {
     125             :           heat2 = 0.0;
     126             :         }
     127       65940 :         _estimate_power(i_pin) += _ref_qprime(i_pin) * (heat1 + heat2) * dz / 2.0;
     128             :       }
     129             :     }
     130             :   }
     131             : 
     132             :   // if a Pin has zero power (_ref_qprime(i_pin) = 0) then I need to avoid dividing by zero. I
     133             :   // divide by a wrong non-zero number which is not correct but this error doesn't mess things cause
     134             :   // _ref_qprime(i_pin) = 0.0
     135         134 :   auto total_power = 0.0;
     136        4744 :   for (unsigned int i_pin = 0; i_pin < n_pins; i_pin++)
     137             :   {
     138        4610 :     total_power += _estimate_power(i_pin);
     139        4610 :     if (_estimate_power(i_pin) == 0.0)
     140         516 :       _estimate_power(i_pin) = 1.0;
     141             :   }
     142             :   // We need to correct the linear power assigned to the nodes of each pin
     143             :   // so that the total power calculated  by the trapezoidal rule agrees with the power assigned by
     144             :   // the user.
     145         134 :   _pin_power_correction = _ref_power.cwiseQuotient(_estimate_power);
     146         134 :   _console << "###########################################" << std::endl;
     147         134 :   _console << "Total power estimation by IC kernel before correction: " << total_power << " [W] "
     148         134 :            << std::endl;
     149         134 :   _console << "IC Power correction vector :\n" << _pin_power_correction << " \n";
     150         134 : }
     151             : 
     152             : Real
     153      258876 : SCMTriPowerIC::value(const Point & p)
     154             : {
     155             :   auto heat_rate = 0.0;
     156      258876 :   auto heated_length = _mesh.getHeatedLength();
     157      258876 :   auto unheated_length_entry = _mesh.getHeatedLengthEntry();
     158             :   Point p1(0, 0, unheated_length_entry);
     159             :   Point P = p - p1;
     160      258876 :   auto pin_mesh_exist = _mesh.pinMeshExist();
     161             : 
     162             :   /// assign power to the nodes located inside the heated section
     163      258876 :   if (MooseUtils::absoluteFuzzyGreaterEqual(p(2), unheated_length_entry) &&
     164      212382 :       MooseUtils::absoluteFuzzyLessEqual(p(2), unheated_length_entry + heated_length))
     165             :   {
     166      167766 :     if (pin_mesh_exist)
     167             :     {
     168             :       // project axial heat rate on pins
     169       52278 :       auto i_pin = _mesh.getPinIndexFromPoint(p);
     170       52278 :       return _ref_qprime(i_pin) * _pin_power_correction(i_pin) * _axial_heat_rate.value(_t, P);
     171             :     }
     172             :     else
     173             :     {
     174             :       // Determine which subchannel this point is in.
     175      115488 :       auto i_ch = _mesh.getSubchannelIndexFromPoint(p);
     176      115488 :       auto subch_type = _mesh.getSubchannelType(i_ch);
     177             :       // project axial heat rate on subchannels
     178             :       {
     179             :         double factor;
     180      115488 :         switch (subch_type)
     181             :         {
     182             :           case EChannelType::CENTER:
     183             :             factor = 1.0 / 6.0;
     184             :             break;
     185       31968 :           case EChannelType::EDGE:
     186             :             factor = 1.0 / 4.0;
     187       31968 :             break;
     188             :           case EChannelType::CORNER:
     189             :             factor = 1.0 / 6.0;
     190             :             break;
     191             :           default:
     192             :             return 0.0; // handle invalid subch_type if needed
     193             :         }
     194      399456 :         for (auto i_pin : _mesh.getChannelPins(i_ch))
     195             :         {
     196      283968 :           heat_rate += factor * _ref_qprime(i_pin) * _pin_power_correction(i_pin) *
     197      283968 :                        _axial_heat_rate.value(_t, P);
     198             :         }
     199      115488 :         return heat_rate;
     200             :       }
     201             :     }
     202             :   }
     203             :   else
     204             :     return 0.0;
     205             : }