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