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             : 
      17             : InputParameters
      18          42 : SCMQuadPowerAux::validParams()
      19             : {
      20          42 :   InputParameters params = AuxKernel::validParams();
      21          42 :   params.addClassDescription(
      22             :       "Computes axial power rate (W/m) assigned to the fuel pins in a quadrilateral lattice "
      23             :       "arrangement");
      24          84 :   params.addRequiredParam<PostprocessorName>(
      25             :       "power", "The postprocessor or Real to use for the total power of the subassembly [W]");
      26          84 :   params.addRequiredParam<std::string>(
      27             :       "filename", "name of radial power profile .txt file (should be a single column) [UnitLess].");
      28          84 :   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          42 :   return params;
      33           0 : }
      34             : 
      35          22 : SCMQuadPowerAux::SCMQuadPowerAux(const InputParameters & parameters)
      36             :   : AuxKernel(parameters),
      37          44 :     _quadMesh(SCM::getConstMesh<QuadSubChannelMesh>(_mesh)),
      38          22 :     _power(getPostprocessorValue("power")),
      39          22 :     _numberoflines(0),
      40          66 :     _filename(getParam<std::string>("filename")),
      41          44 :     _axial_heat_rate(getFunction("axial_heat_rate"))
      42             : {
      43          22 :   if (processor_id() > 0)
      44           6 :     return;
      45             : 
      46          16 :   if (!_quadMesh.pinMeshExist())
      47           0 :     mooseError(name(), ": This object requires a pin mesh.");
      48             : 
      49          16 :   auto n_pins = _quadMesh.getNumOfPins();
      50             :   // Matrix sizing
      51          16 :   _power_dis.resize(n_pins, 1);
      52          16 :   _power_dis.setZero();
      53          16 :   _pin_power_correction.resize(n_pins, 1);
      54          16 :   _pin_power_correction.setOnes();
      55             : 
      56             :   Real vin;
      57          16 :   std::ifstream inFile;
      58             : 
      59          16 :   inFile.open(_filename);
      60          16 :   if (!inFile)
      61           0 :     mooseError(name(), "unable to open file : ", _filename);
      62             : 
      63         164 :   while (inFile >> vin)
      64         148 :     _numberoflines += 1;
      65             : 
      66          16 :   if (inFile.fail() && !inFile.eof())
      67           0 :     mooseError(name(), " non numerical input at line : ", _numberoflines);
      68             : 
      69          16 :   if (_numberoflines != n_pins)
      70           0 :     mooseError(name(), " Radial profile file doesn't have correct size : ", n_pins);
      71          16 :   inFile.close();
      72             : 
      73          16 :   inFile.open(_filename);
      74             :   int i(0);
      75         164 :   while (inFile >> vin)
      76             :   {
      77         148 :     _power_dis(i, 0) = vin;
      78         148 :     i++;
      79             :   }
      80          16 :   inFile.close();
      81          16 :   _console << " Power distribution matrix :\n" << _power_dis << " \n";
      82          16 : }
      83             : 
      84             : void
      85          44 : SCMQuadPowerAux::initialSetup()
      86             : {
      87          44 :   if (processor_id() > 0)
      88          12 :     return;
      89             : 
      90          32 :   auto n_pins = _quadMesh.getNumOfPins();
      91          32 :   auto nz = _quadMesh.getNumOfAxialCells();
      92          32 :   auto z_grid = _quadMesh.getZGrid();
      93          32 :   auto heated_length = _quadMesh.getHeatedLength();
      94          32 :   auto unheated_length_entry = _quadMesh.getHeatedLengthEntry();
      95          32 :   auto sum = _power_dis.sum();
      96             : 
      97             :   // full (100%) power of one pin [W]
      98          32 :   auto fpin_power = _power / sum;
      99             :   // actual pin power [W]
     100          32 :   _ref_power = _power_dis * fpin_power;
     101             :   // Convert the actual pin power to a linear heat rate [W/m]
     102          32 :   _ref_qprime = _ref_power / heated_length;
     103             : 
     104          32 :   _estimate_power.resize(n_pins, 1);
     105          32 :   _estimate_power.setZero();
     106         424 :   for (unsigned int iz = 1; iz < nz + 1; iz++)
     107             :   {
     108             :     // Compute axial location of nodes.
     109         392 :     auto z2 = z_grid[iz];
     110         392 :     auto z1 = z_grid[iz - 1];
     111         392 :     Point p1(0, 0, z1 - unheated_length_entry);
     112         392 :     Point p2(0, 0, z2 - unheated_length_entry);
     113         392 :     auto heat1 = _axial_heat_rate.value(_t, p1);
     114         392 :     auto heat2 = _axial_heat_rate.value(_t, p2);
     115         392 :     if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
     116         296 :         MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
     117             :     {
     118             :       // cycle through pins
     119        2680 :       for (unsigned int i_pin = 0; i_pin < n_pins; i_pin++)
     120             :       {
     121             :         // Compute the height of this element.
     122        2480 :         auto dz = z2 - z1;
     123             : 
     124             :         // calculation of power for the first heated segment if nodes don't align
     125        2480 :         if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
     126             :             MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry))
     127             :         {
     128             :           heat1 = 0.0;
     129             :         }
     130             : 
     131             :         // calculation of power for the last heated segment if nodes don't align
     132        2480 :         if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry + heated_length) &&
     133             :             MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
     134             :         {
     135             :           heat2 = 0.0;
     136             :         }
     137             : 
     138        2480 :         _estimate_power(i_pin) += _ref_qprime(i_pin) * (heat1 + heat2) * dz / 2.0;
     139             :       }
     140             :     }
     141             :   }
     142             : 
     143             :   // if a Pin has zero power (_ref_qprime(i_pin) = 0) then I need to avoid dividing by zero. I
     144             :   // divide by a wrong non-zero number which is not correct but this error doesn't mess things
     145             :   // cause _ref_qprime(i_pin) = 0.0
     146          32 :   auto total_power = 0.0;
     147         328 :   for (unsigned int i_pin = 0; i_pin < n_pins; i_pin++)
     148             :   {
     149         296 :     total_power += _estimate_power(i_pin);
     150         296 :     if (_estimate_power(i_pin) == 0.0)
     151           0 :       _estimate_power(i_pin) = 1.0;
     152             :   }
     153             :   // We need to correct the linear power assigned to the nodes of each pin
     154             :   // so that the total power calculated  by the trapezoidal rule agrees with the power assigned
     155             :   // by the user.
     156          32 :   _pin_power_correction = _ref_power.cwiseQuotient(_estimate_power);
     157          32 :   _console << "###########################################" << std::endl;
     158          32 :   _console << "Total power estimation by Aux kernel before correction: " << total_power << " [W] "
     159          32 :            << std::endl;
     160          32 :   _console << "Aux Power correction vector :\n" << _pin_power_correction << " \n";
     161          32 : }
     162             : 
     163             : Real
     164        2760 : SCMQuadPowerAux::computeValue()
     165             : {
     166        2760 :   Point p = *_current_node;
     167        2760 :   auto heated_length = _quadMesh.getHeatedLength();
     168        2760 :   auto unheated_length_entry = _quadMesh.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        2760 :   if (MooseUtils::absoluteFuzzyGreaterEqual(p(2), unheated_length_entry) &&
     174        2580 :       MooseUtils::absoluteFuzzyLessEqual(p(2), unheated_length_entry + heated_length))
     175             :   {
     176        2400 :     auto i_pin = _quadMesh.getPinIndexFromPoint(p);
     177        2400 :     return _ref_qprime(i_pin) * _pin_power_correction(i_pin) * _axial_heat_rate.value(_t, P);
     178             :   }
     179             :   else
     180             :     return 0.0;
     181             : }