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 "QuadInterWrapperPowerIC.h"
      11             : #include "Function.h"
      12             : #include "QuadInterWrapperMesh.h"
      13             : 
      14             : using namespace std;
      15             : using namespace Eigen;
      16             : 
      17             : registerMooseObject("SubChannelApp", QuadInterWrapperPowerIC);
      18             : 
      19             : InputParameters
      20          68 : QuadInterWrapperPowerIC::validParams()
      21             : {
      22          68 :   InputParameters params = QuadInterWrapperBaseIC::validParams();
      23          68 :   params.addClassDescription("Computes axial power rate, W/m that goes into the inter-wrapper "
      24             :                              "cells in a square lattice subchannel arrangement");
      25         136 :   params.addParam<Real>(
      26         136 :       "power", 0.0, "The power of all the sub-assemblies that the inter-wrapper wraps around[W]");
      27         136 :   params.addParam<std::string>("filename",
      28             :                                "file_was_not_found",
      29             :                                "name of power profile .txt file (should be a single column). It's "
      30             :                                "a Radial Power Profile. [UnitLess]");
      31         136 :   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             :       "Zero is considered the inlet of the heated length.");
      37          68 :   return params;
      38           0 : }
      39             : 
      40          35 : QuadInterWrapperPowerIC::QuadInterWrapperPowerIC(const InputParameters & params)
      41             :   : QuadInterWrapperBaseIC(params),
      42          35 :     _power(getParam<Real>("power")),
      43          35 :     _numberoflines(0),
      44         105 :     _filename(getParam<std::string>("filename")),
      45          70 :     _axial_heat_rate(getFunction("axial_heat_rate"))
      46             : {
      47          35 :   auto nx = _mesh.getNx();
      48          35 :   auto ny = _mesh.getNy();
      49          35 :   auto heated_length = _mesh.getHeatedLength();
      50             : 
      51          35 :   _power_dis.resize((ny - 1) * (nx - 1), 1);
      52          35 :   _power_dis.setZero();
      53          35 :   _assembly_power_correction.resize((ny - 1) * (nx - 1), 1);
      54          35 :   _assembly_power_correction.setOnes();
      55             :   Real vin;
      56          35 :   ifstream inFile;
      57             : 
      58          35 :   _console << "Power file: " << _filename << std::endl;
      59             : 
      60          35 :   if (_filename.compare("file_was_not_found"))
      61             :   {
      62           0 :     inFile.open(_filename);
      63           0 :     if (!inFile)
      64           0 :       mooseError(name(), "unable to open file : ", _filename);
      65             : 
      66           0 :     while (inFile >> vin)
      67           0 :       _numberoflines += 1;
      68             : 
      69           0 :     if (inFile.fail() && !inFile.eof())
      70           0 :       mooseError(name(), " non numerical input at line : ", _numberoflines);
      71             : 
      72           0 :     if (_numberoflines != (ny - 1) * (nx - 1))
      73           0 :       mooseError(name(), " Radial profile file doesn't have correct size : ", (ny - 1) * (nx - 1));
      74           0 :     inFile.close();
      75             :   }
      76             :   else
      77             :   {
      78          35 :     _numberoflines = (ny - 1) * (nx - 1);
      79             :   }
      80             : 
      81          35 :   if (_filename.compare("file_was_not_found"))
      82             :   {
      83           0 :     inFile.open(_filename);
      84             :     int i(0);
      85           0 :     while (inFile >> vin)
      86             :     {
      87           0 :       _power_dis(i, 0) = vin;
      88           0 :       i++;
      89             :     }
      90           0 :     inFile.close();
      91             :   }
      92             : 
      93          35 :   _console << " Power distribution matrix :\n" << _power_dis << " \n";
      94          35 :   auto sum = _power_dis.sum();
      95             :   // full (100%) power of one pin [W]
      96          35 :   auto fpin_power = _power / sum;
      97             :   // actual pin power [W]
      98          35 :   _ref_power = _power_dis * fpin_power;
      99             :   // Convert the actual pin power to a linear heat rate [W/m]
     100          35 :   _ref_qprime = _ref_power / heated_length;
     101          35 : }
     102             : 
     103             : void
     104          24 : QuadInterWrapperPowerIC::initialSetup()
     105             : {
     106             : 
     107          24 :   auto nx = _mesh.getNx();
     108          24 :   auto ny = _mesh.getNy();
     109          24 :   auto nz = _mesh.getNumOfAxialCells();
     110          24 :   auto z_grid = _mesh.getZGrid();
     111          24 :   auto heated_length = _mesh.getHeatedLength();
     112          24 :   auto unheated_length_entry = _mesh.getHeatedLengthEntry();
     113             : 
     114          24 :   _estimate_power.resize((ny - 1) * (nx - 1), 1);
     115          24 :   _estimate_power.setZero();
     116         264 :   for (unsigned int iz = 1; iz < nz + 1; iz++)
     117             :   {
     118             :     // Compute the height of this element.
     119         240 :     auto dz = z_grid[iz] - z_grid[iz - 1];
     120             :     // Compute axial location of nodes.
     121             :     auto z2 = z_grid[iz];
     122             :     auto z1 = z_grid[iz - 1];
     123         240 :     Point p1(0, 0, z1 - unheated_length_entry);
     124         240 :     Point p2(0, 0, z2 - unheated_length_entry);
     125             :     // cycle through pins to estimate the total power of each pin
     126         240 :     if (z2 > unheated_length_entry && z2 <= unheated_length_entry + heated_length)
     127             :     {
     128        6240 :       for (unsigned int i_pin = 0; i_pin < (ny - 1) * (nx - 1); i_pin++)
     129             :       {
     130             :         // use of trapezoidal rule to add to total power of pin
     131        6000 :         _estimate_power(i_pin) +=
     132        6000 :             _ref_qprime(i_pin) * (_axial_heat_rate.value(_t, p1) + _axial_heat_rate.value(_t, p2)) *
     133        6000 :             dz / 2.0;
     134             :       }
     135             :     }
     136             :   }
     137             : 
     138             :   // if a Pin has zero power (_ref_qprime(j, i) = 0) then I need to avoid dividing by zero. I
     139             :   // divide by a wrong non-zero number which is not correct but this error doesn't mess things cause
     140             :   // _ref_qprime(j, i) = 0.0
     141         624 :   for (unsigned int i_pin = 0; i_pin < (ny - 1) * (nx - 1); i_pin++)
     142             :   {
     143         600 :     if (_estimate_power(i_pin) == 0.0)
     144           0 :       _estimate_power(i_pin) = 1.0;
     145             :   }
     146          24 :   _assembly_power_correction = _ref_power.cwiseQuotient(_estimate_power);
     147          24 : }
     148             : 
     149             : Real
     150       17280 : QuadInterWrapperPowerIC::value(const Point & p)
     151             : {
     152       17280 :   auto heated_length = _mesh.getHeatedLength();
     153       17280 :   auto unheated_length_entry = _mesh.getHeatedLengthEntry();
     154             :   Point p1(0, 0, unheated_length_entry);
     155             :   Point P = p - p1;
     156       17280 :   auto pin_mesh_exist = _mesh.pinMeshExist();
     157             : 
     158       17280 :   if (pin_mesh_exist)
     159             :   {
     160             :     // project axial heat rate on pins
     161           0 :     auto i_pin = _mesh.getPinIndexFromPoint(p);
     162             :     {
     163           0 :       if (p(2) >= unheated_length_entry && p(2) <= unheated_length_entry + heated_length)
     164           0 :         return _ref_qprime(i_pin) * _assembly_power_correction(i_pin) *
     165           0 :                _axial_heat_rate.value(_t, P);
     166             :       else
     167             :         return 0.0;
     168             :     }
     169             :   }
     170             :   else
     171             :   {
     172             :     // project axial heat rate on subchannels
     173       17280 :     auto i_ch = _mesh.getSubchannelIndexFromPoint(p);
     174             :     // if we are adjacent to the heated part of the fuel Pin
     175       17280 :     if (p(2) >= unheated_length_entry && p(2) <= unheated_length_entry + heated_length)
     176             :     {
     177             :       auto heat_rate = 0.0;
     178       65280 :       for (auto i_pin : _mesh.getChannelPins(i_ch))
     179             :       {
     180       48000 :         heat_rate += 0.25 * _ref_qprime(i_pin) * _assembly_power_correction(i_pin) *
     181       48000 :                      _axial_heat_rate.value(_t, P);
     182             :       }
     183       17280 :       return heat_rate;
     184             :     }
     185             :     else
     186             :       return 0.0;
     187             :   }
     188             : }