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