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             : 
      19             : InputParameters
      20         238 : SCMQuadPowerIC::validParams()
      21             : {
      22         238 :   InputParameters params = QuadSubChannelBaseIC::validParams();
      23         238 :   params.addClassDescription("Computes axial heat rate (W/m) assigned to the fuel pins in a square "
      24             :                              "lattice arrangement");
      25             :   // params.addRequiredParam<Real>("power", "The total power of the subassembly [W]");
      26         476 :   params.addRequiredParam<PostprocessorName>(
      27             :       "power", "The postprocessor or Real to use for the total power of the subassembly [W]");
      28         476 :   params.addRequiredParam<std::string>(
      29             :       "filename", "name of radial power profile .txt file (should be a single column) [UnitLess].");
      30         476 :   params.addParam<FunctionName>(
      31             :       "axial_heat_rate",
      32             :       "1.0",
      33             :       "user provided normalized function of axial heat rate [Unitless]. "
      34             :       "The integral over pin heated length should equal the heated length");
      35         238 :   return params;
      36           0 : }
      37             : 
      38         123 : SCMQuadPowerIC::SCMQuadPowerIC(const InputParameters & params)
      39             :   : QuadSubChannelBaseIC(params),
      40         123 :     _power(getPostprocessorValue("power")),
      41         123 :     _numberoflines(0),
      42         369 :     _filename(getParam<std::string>("filename")),
      43         246 :     _axial_heat_rate(getFunction("axial_heat_rate"))
      44             : {
      45         123 :   if (processor_id() > 0)
      46          33 :     return;
      47             : 
      48          90 :   if (!_mesh.pinMeshExist())
      49           0 :     mooseError(name(), ": This object requires a pin mesh.");
      50             : 
      51          90 :   auto n_pins = _mesh.getNumOfPins();
      52          90 :   auto heated_length = _mesh.getHeatedLength();
      53             : 
      54          90 :   _power_dis.resize(n_pins, 1);
      55          90 :   _power_dis.setZero();
      56          90 :   _pin_power_correction.resize(n_pins, 1);
      57          90 :   _pin_power_correction.setOnes();
      58             : 
      59             :   Real vin;
      60          90 :   ifstream inFile;
      61             : 
      62          90 :   inFile.open(_filename);
      63          90 :   if (!inFile)
      64           0 :     mooseError(name(), "unable to open file : ", _filename);
      65             : 
      66        1138 :   while (inFile >> vin)
      67        1048 :     _numberoflines += 1;
      68             : 
      69          90 :   if (inFile.fail() && !inFile.eof())
      70           0 :     mooseError(name(), " non numerical input at line : ", _numberoflines);
      71             : 
      72          90 :   if (_numberoflines != n_pins)
      73           0 :     mooseError(name(), " Radial profile file doesn't have correct size : ", n_pins);
      74          90 :   inFile.close();
      75             : 
      76          90 :   inFile.open(_filename);
      77             :   int i(0);
      78        1138 :   while (inFile >> vin)
      79             :   {
      80        1048 :     _power_dis(i, 0) = vin;
      81        1048 :     i++;
      82             :   }
      83          90 :   inFile.close();
      84          90 :   _console << " Power distribution matrix :\n" << _power_dis << " \n";
      85             : 
      86          90 :   auto sum = _power_dis.sum();
      87             :   // full (100%) power of one pin [W]
      88          90 :   auto fpin_power = _power / sum;
      89             :   // actual pin power [W]
      90          90 :   _ref_power = _power_dis * fpin_power;
      91             :   // Convert the actual pin power to a linear heat rate [W/m]
      92          90 :   _ref_qprime = _ref_power / heated_length;
      93          90 : }
      94             : 
      95             : void
      96         108 : SCMQuadPowerIC::initialSetup()
      97             : {
      98         108 :   if (processor_id() > 0)
      99          31 :     return;
     100          77 :   auto n_pins = _mesh.getNumOfPins();
     101          77 :   auto nz = _mesh.getNumOfAxialCells();
     102          77 :   auto z_grid = _mesh.getZGrid();
     103          77 :   auto heated_length = _mesh.getHeatedLength();
     104          77 :   auto unheated_length_entry = _mesh.getHeatedLengthEntry();
     105             : 
     106          77 :   _estimate_power.resize(n_pins, 1);
     107          77 :   _estimate_power.setZero();
     108        1485 :   for (unsigned int iz = 1; iz < nz + 1; iz++)
     109             :   {
     110             :     // Compute axial location of nodes.
     111        1408 :     auto z2 = z_grid[iz];
     112        1408 :     auto z1 = z_grid[iz - 1];
     113        1408 :     Point p1(0, 0, z1 - unheated_length_entry);
     114        1408 :     Point p2(0, 0, z2 - unheated_length_entry);
     115        1408 :     auto heat1 = _axial_heat_rate.value(_t, p1);
     116        1408 :     auto heat2 = _axial_heat_rate.value(_t, p2);
     117        1408 :     if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
     118        1297 :         MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
     119             :     {
     120             :       // cycle through pins
     121       13052 :       for (unsigned int i_pin = 0; i_pin < n_pins; i_pin++)
     122             :       {
     123             :         // Compute the height of this element.
     124       11818 :         auto dz = z2 - z1;
     125             : 
     126             :         // calculation of power for the first heated segment if nodes don't align
     127       11818 :         if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
     128             :             MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry))
     129             :         {
     130             :           heat1 = 0.0;
     131             :         }
     132             : 
     133             :         // calculation of power for the last heated segment if nodes don't align
     134       11818 :         if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry + heated_length) &&
     135             :             MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
     136             :         {
     137             :           heat2 = 0.0;
     138             :         }
     139             : 
     140       11818 :         _estimate_power(i_pin) += _ref_qprime(i_pin) * (heat1 + heat2) * dz / 2.0;
     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 cause
     146             :   // _ref_qprime(i_pin) = 0.0
     147          77 :   auto total_power = 0.0;
     148         800 :   for (unsigned int i_pin = 0; i_pin < n_pins; i_pin++)
     149             :   {
     150         723 :     total_power += _estimate_power(i_pin);
     151         723 :     if (_estimate_power(i_pin) == 0.0)
     152          72 :       _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 by
     156             :   // the user.
     157          77 :   _pin_power_correction = _ref_power.cwiseQuotient(_estimate_power);
     158          77 :   _console << "###########################################" << std::endl;
     159          77 :   _console << "Total power estimation by IC kernel before correction: " << total_power << " [W] "
     160          77 :            << std::endl;
     161          77 :   _console << "IC Power correction vector :\n" << _pin_power_correction << " \n";
     162          77 : }
     163             : 
     164             : Real
     165       25448 : SCMQuadPowerIC::value(const Point & p)
     166             : {
     167       25448 :   auto heated_length = _mesh.getHeatedLength();
     168       25448 :   auto unheated_length_entry = _mesh.getHeatedLengthEntry();
     169             :   Point p1(0, 0, unheated_length_entry);
     170             :   Point P = p - p1;
     171             : 
     172             :   /// assign power to the nodes located within the heated section
     173       25448 :   if (MooseUtils::absoluteFuzzyGreaterEqual(p(2), unheated_length_entry) &&
     174       23916 :       MooseUtils::absoluteFuzzyLessEqual(p(2), unheated_length_entry + heated_length))
     175             :   {
     176       23488 :     auto i_pin = _mesh.getPinIndexFromPoint(p);
     177       23488 :     return _ref_qprime(i_pin) * _pin_power_correction(i_pin) * _axial_heat_rate.value(_t, P);
     178             :   }
     179             :   else
     180             :     return 0.0;
     181             : }