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