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 "TriInterWrapperPowerIC.h"
      11             : #include "Function.h"
      12             : #include "TriInterWrapperMesh.h"
      13             : 
      14             : registerMooseObject("SubChannelApp", TriInterWrapperPowerIC);
      15             : 
      16             : InputParameters
      17          80 : TriInterWrapperPowerIC::validParams()
      18             : {
      19          80 :   InputParameters params = TriInterWrapperBaseIC::validParams();
      20          80 :   params.addClassDescription("Computes linear power rate (W/m) that goes into interwrapper cells "
      21             :                              "in a triangular subchannel lattice");
      22         160 :   params.addParam<Real>("power", 0.0, "The total heating power [W]");
      23         160 :   params.addParam<std::string>("filename",
      24             :                                "file_was_not_found",
      25             :                                "name of power profile .txt file (should be a single column). It's "
      26             :                                "a Radial Power Profile. [UnitLess]");
      27         160 :   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          80 :   return params;
      32           0 : }
      33             : 
      34          41 : TriInterWrapperPowerIC::TriInterWrapperPowerIC(const InputParameters & params)
      35             :   : TriInterWrapperBaseIC(params),
      36          41 :     _power(getParam<Real>("power")),
      37          41 :     _numberoflines(0),
      38         123 :     _filename(getParam<std::string>("filename")),
      39          82 :     _axial_heat_rate(getFunction("axial_heat_rate"))
      40             : {
      41          41 :   auto n_assemblies = _mesh.getNumOfAssemblies();
      42             : 
      43          41 :   _power_dis.resize(n_assemblies);
      44          41 :   _pin_power_correction.resize(n_assemblies);
      45          41 :   _ref_power.resize(n_assemblies);
      46          41 :   _ref_qprime.resize(n_assemblies);
      47          41 :   _estimate_power.resize(n_assemblies);
      48        1018 :   for (unsigned int i = 0; i < n_assemblies; i++)
      49             :   {
      50         977 :     _power_dis[i] = 0.0;
      51         977 :     _pin_power_correction[i] = 1.0;
      52             :   }
      53             : 
      54          41 :   if (_filename.compare("file_was_not_found"))
      55             :   {
      56             : 
      57             :     Real vin;
      58           0 :     std::ifstream inFile;
      59             : 
      60           0 :     while (inFile >> vin)
      61           0 :       _numberoflines += 1;
      62             : 
      63           0 :     inFile.open(_filename);
      64           0 :     if (!inFile)
      65           0 :       mooseError(name(), ": Unable to open file: ", _filename);
      66             : 
      67             :     if (inFile.fail() && !inFile.eof())
      68             :       mooseError(name(), ": Non numerical input at line: ", _numberoflines);
      69             : 
      70           0 :     if (_numberoflines != n_assemblies)
      71           0 :       mooseError(name(), ": Radial profile file doesn't have correct size: ", n_assemblies);
      72           0 :     inFile.close();
      73             : 
      74           0 :     inFile.open(_filename);
      75             :     int i = 0;
      76           0 :     while (inFile >> vin)
      77             :     {
      78           0 :       _power_dis[i] = vin;
      79           0 :       i++;
      80             :     }
      81           0 :     inFile.close();
      82           0 :   }
      83             :   Real sum = 0.0;
      84             : 
      85        1018 :   for (unsigned int i = 0; i < n_assemblies; i++)
      86             :   {
      87         977 :     sum = sum + _power_dis[i];
      88             :   }
      89             :   // full pin (100%) power of one pin [W]
      90          41 :   auto fpin_power = _power / sum;
      91             :   // actual pin power [W]
      92        1018 :   for (unsigned int i = 0; i < n_assemblies; i++)
      93             :   {
      94         977 :     _ref_power[i] = _power_dis[i] * fpin_power;
      95             :   }
      96             : 
      97             :   // Convert the actual pin power to a linear heat rate [W/m]
      98          41 :   auto heated_length = _mesh.getHeatedLength();
      99             : 
     100        1018 :   for (unsigned int i = 0; i < n_assemblies; i++)
     101             :   {
     102         977 :     _ref_qprime[i] = _ref_power[i] / heated_length;
     103             :   }
     104          41 : }
     105             : 
     106             : void
     107          30 : TriInterWrapperPowerIC::initialSetup()
     108             : {
     109          30 :   auto n_assemblies = _mesh.getNumOfAssemblies();
     110          30 :   auto nz = _mesh.getNumOfAxialCells();
     111          30 :   auto z_grid = _mesh.getZGrid();
     112          30 :   auto heated_length = _mesh.getHeatedLength();
     113          30 :   auto unheated_length_entry = _mesh.getHeatedLengthEntry();
     114             : 
     115          30 :   _estimate_power.resize(n_assemblies);
     116             : 
     117         330 :   for (unsigned int iz = 1; iz < nz + 1; iz++)
     118             :   {
     119             :     // Compute the height of this element.
     120         300 :     auto dz = z_grid[iz] - z_grid[iz - 1];
     121             :     // Compute axial location of nodes.
     122             :     auto z2 = z_grid[iz];
     123             :     auto z1 = z_grid[iz - 1];
     124         300 :     Point p1(0, 0, z1 - unheated_length_entry);
     125         300 :     Point p2(0, 0, z2 - unheated_length_entry);
     126             :     // cycle through pins
     127             : 
     128         300 :     if (z2 > unheated_length_entry && z2 <= unheated_length_entry + heated_length)
     129             :     {
     130        6000 :       for (unsigned int i_pin = 0; i_pin < n_assemblies; i_pin++)
     131             :       {
     132             :         // use of trapezoidal rule  to calculate local power
     133        5700 :         _estimate_power[i_pin] +=
     134        5700 :             _ref_qprime[i_pin] * (_axial_heat_rate.value(_t, p1) + _axial_heat_rate.value(_t, p2)) *
     135        5700 :             dz / 2.0;
     136             :       }
     137             :     }
     138             :   }
     139         600 :   for (unsigned int i_pin = 0; i_pin < n_assemblies; i_pin++)
     140             :   {
     141         570 :     _pin_power_correction[i_pin] = _ref_power[i_pin] / _estimate_power[i_pin];
     142             :   }
     143          30 : }
     144             : 
     145             : Real
     146       25200 : TriInterWrapperPowerIC::value(const Point & p)
     147             : {
     148             :   // Determine which subchannel this point is in.
     149       25200 :   auto i = _mesh.getSubchannelIndexFromPoint(p);
     150       25200 :   auto subch_type = _mesh.getSubchannelType(i);
     151             :   Real sum = 0.0;
     152       25200 :   auto heated_length = _mesh.getHeatedLength();
     153       25200 :   auto unheated_length_entry = _mesh.getHeatedLengthEntry();
     154             :   Point p1(0, 0, unheated_length_entry);
     155             :   // Point P = p - p1;
     156             : 
     157             :   // if we are adjacent to the heated part of the fuel Pin
     158       25200 :   if (p(2) >= unheated_length_entry && p(2) <= unheated_length_entry + heated_length)
     159             :   {
     160       25200 :     if (subch_type == EChannelType::CENTER)
     161             :     {
     162       57600 :       for (unsigned int j = 0; j < 3; j++)
     163             :       {
     164       43200 :         auto rod_idx = _mesh.getPinIndex(i, j);
     165       43200 :         sum = sum + 1.0 / 6.0 * _ref_qprime[rod_idx] * _pin_power_correction[rod_idx] *
     166       43200 :                         _axial_heat_rate.value(_t, p);
     167             :       }
     168       14400 :       return sum;
     169             :     }
     170       10800 :     else if (subch_type == EChannelType::EDGE)
     171             :     {
     172       21600 :       for (unsigned int j = 0; j < 2; j++)
     173             :       {
     174       14400 :         auto rod_idx = _mesh.getPinIndex(i, j);
     175       14400 :         sum = sum + 1.0 / 4.0 * _ref_qprime[rod_idx] * _pin_power_correction[rod_idx] *
     176       14400 :                         _axial_heat_rate.value(_t, p);
     177             :       }
     178        7200 :       return sum;
     179             :     }
     180        3600 :     else if (subch_type == EChannelType::CORNER)
     181             :     {
     182        3600 :       auto rod_idx = _mesh.getPinIndex(i, 0);
     183        3600 :       sum = 1.0 / 6.0 * _ref_qprime[rod_idx] * _pin_power_correction[rod_idx] *
     184        3600 :             _axial_heat_rate.value(_t, p);
     185        3600 :       return sum;
     186             :     }
     187             :   }
     188             :   return 0.;
     189             : }