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