https://mooseframework.inl.gov
SCMQuadPowerIC.C
Go to the documentation of this file.
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 registerMooseObjectRenamed("SubChannelApp", QuadPowerIC, "06/30/2025 24:00", SCMQuadPowerIC);
19 
22 {
24  params.addClassDescription("Computes axial heat rate (W/m) that goes into the subchannel cells "
25  "or is assigned to the fuel pins, in a square lattice arrangement");
26  // params.addRequiredParam<Real>("power", "The total power of the subassembly [W]");
27  params.addRequiredParam<PostprocessorName>(
28  "power", "The postprocessor or Real to use for the total power of the subassembly [W]");
29  params.addRequiredParam<std::string>(
30  "filename", "name of radial power profile .txt file (should be a single column) [UnitLess].");
31  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  return params;
37 }
38 
40  : QuadSubChannelBaseIC(params),
41  _power(getPostprocessorValue("power")),
42  _numberoflines(0),
43  _filename(getParam<std::string>("filename")),
44  _axial_heat_rate(getFunction("axial_heat_rate"))
45 {
46  auto nx = _mesh.getNx();
47  auto ny = _mesh.getNy();
48  auto heated_length = _mesh.getHeatedLength();
49 
50  _power_dis.resize((ny - 1) * (nx - 1), 1);
51  _power_dis.setZero();
52  _pin_power_correction.resize((ny - 1) * (nx - 1), 1);
53  _pin_power_correction.setOnes();
54 
55  Real vin;
56  ifstream inFile;
57 
58  inFile.open(_filename);
59  if (!inFile)
60  mooseError(name(), "unable to open file : ", _filename);
61 
62  while (inFile >> vin)
63  _numberoflines += 1;
64 
65  if (inFile.fail() && !inFile.eof())
66  mooseError(name(), " non numerical input at line : ", _numberoflines);
67 
68  if (_numberoflines != (ny - 1) * (nx - 1))
69  mooseError(name(), " Radial profile file doesn't have correct size : ", (ny - 1) * (nx - 1));
70  inFile.close();
71 
72  inFile.open(_filename);
73  int i(0);
74  while (inFile >> vin)
75  {
76  _power_dis(i, 0) = vin;
77  i++;
78  }
79  inFile.close();
80  _console << " Power distribution matrix :\n" << _power_dis << " \n";
81 
82  auto sum = _power_dis.sum();
83  // full (100%) power of one pin [W]
84  auto fpin_power = _power / sum;
85  // actual pin power [W]
86  _ref_power = _power_dis * fpin_power;
87  // Convert the actual pin power to a linear heat rate [W/m]
88  _ref_qprime = _ref_power / heated_length;
89 }
90 
91 void
93 {
94  auto nx = _mesh.getNx();
95  auto ny = _mesh.getNy();
96  auto n_pins = (nx - 1) * (ny - 1);
97  auto nz = _mesh.getNumOfAxialCells();
98  auto z_grid = _mesh.getZGrid();
99  auto heated_length = _mesh.getHeatedLength();
100  auto unheated_length_entry = _mesh.getHeatedLengthEntry();
101 
102  _estimate_power.resize(n_pins, 1);
103  _estimate_power.setZero();
104  for (unsigned int iz = 1; iz < nz + 1; iz++)
105  {
106  // Compute axial location of nodes.
107  auto z2 = z_grid[iz];
108  auto z1 = z_grid[iz - 1];
109  Point p1(0, 0, z1 - unheated_length_entry);
110  Point p2(0, 0, z2 - unheated_length_entry);
111  auto heat1 = _axial_heat_rate.value(_t, p1);
112  auto heat2 = _axial_heat_rate.value(_t, p2);
113  if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
114  MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
115  {
116  // cycle through pins
117  for (unsigned int i_pin = 0; i_pin < n_pins; i_pin++)
118  {
119  // Compute the height of this element.
120  auto dz = z2 - z1;
121 
122  // calculation of power for the first heated segment if nodes don't align
123  if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
124  MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry))
125  {
126  heat1 = 0.0;
127  }
128 
129  // calculation of power for the last heated segment if nodes don't align
130  if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry + heated_length) &&
131  MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
132  {
133  heat2 = 0.0;
134  }
135 
136  _estimate_power(i_pin) += _ref_qprime(i_pin) * (heat1 + heat2) * dz / 2.0;
137  }
138  }
139  }
140  // if a Pin has zero power (_ref_qprime(i_pin) = 0) then I need to avoid dividing by zero. I
141  // divide by a wrong non-zero number which is not correct but this error doesn't mess things cause
142  // _ref_qprime(i_pin) = 0.0
143  auto total_power = 0.0;
144  for (unsigned int i_pin = 0; i_pin < n_pins; i_pin++)
145  {
146  total_power += _estimate_power(i_pin);
147  if (_estimate_power(i_pin) == 0.0)
148  _estimate_power(i_pin) = 1.0;
149  }
150  // We need to correct the linear power assigned to the nodes of each pin
151  // so that the total power calculated by the trapezoidal rule agrees with the power assigned by
152  // the user.
154  _console << "###########################################" << std::endl;
155  _console << "Total power estimation by IC kernel before correction: " << total_power << " [W] "
156  << std::endl;
157  _console << "IC Power correction vector :\n" << _pin_power_correction << " \n";
158 }
159 
160 Real
161 SCMQuadPowerIC::value(const Point & p)
162 {
163  auto heated_length = _mesh.getHeatedLength();
164  auto unheated_length_entry = _mesh.getHeatedLengthEntry();
165  Point p1(0, 0, unheated_length_entry);
166  Point P = p - p1;
167  auto pin_mesh_exist = _mesh.pinMeshExist();
168 
170  if (MooseUtils::absoluteFuzzyGreaterEqual(p(2), unheated_length_entry) &&
171  MooseUtils::absoluteFuzzyLessEqual(p(2), unheated_length_entry + heated_length))
172  {
173  if (pin_mesh_exist)
174  {
175  // project axial heat rate on pins
176  auto i_pin = _mesh.getPinIndexFromPoint(p);
177  return _ref_qprime(i_pin) * _pin_power_correction(i_pin) * _axial_heat_rate.value(_t, P);
178  }
179  else
180  {
181  // project axial heat rate on subchannels
182  auto i_ch = _mesh.getSubchannelIndexFromPoint(p);
183  // if we are adjacent to the heated part of the fuel Pin
184  auto heat_rate = 0.0;
185  for (auto i_pin : _mesh.getChannelPins(i_ch))
186  {
187  heat_rate += 0.25 * _ref_qprime(i_pin) * _pin_power_correction(i_pin) *
189  }
190  return heat_rate;
191  }
192  }
193  else
194  return 0.0;
195 }
SCMQuadPowerIC(const InputParameters &params)
const PostprocessorValue & _power
The total power of the assembly.
registerMooseObject("SubChannelApp", SCMQuadPowerIC)
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
unsigned int _numberoflines
The number of lines associated with the radial power profile .txt file.
virtual const std::vector< Real > & getZGrid() const
Get axial location of layers.
static InputParameters validParams()
An abstract class for ICs for quadrilateral subchannels.
Eigen::MatrixXd _ref_qprime
Average linear heat rate over the whole pin [W/m].
virtual const std::string & name() const
void addRequiredParam(const std::string &name, const std::string &doc_string)
Real value(const Point &p) override
const QuadSubChannelMesh & _mesh
static InputParameters validParams()
unsigned int getSubchannelIndexFromPoint(const Point &p) const override
Return a subchannel index for a given physical point p
virtual const Real & getHeatedLength() const
Return heated length.
Eigen::MatrixXd _estimate_power
Matrix which will hold the total estimated power of each pin [W].
virtual bool pinMeshExist() const override
Return if Pin Mesh exists or not.
Eigen::MatrixXd _pin_power_correction
The correction that will be applied to the estimated calculation [unitless].
bool absoluteFuzzyLessThan(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
Eigen::MatrixXd _power_dis
matrix that holds the values of the relative pin power
Eigen::MatrixXd _ref_power
Actual pin power [W].
virtual const unsigned int & getNumOfAxialCells() const
Return the number of axial cells.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual const unsigned int & getNy() const
Number of subchannels in the -y direction.
bool absoluteFuzzyLessEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
const Function & _axial_heat_rate
bool absoluteFuzzyGreaterEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
unsigned int getPinIndexFromPoint(const Point &p) const override
Return a pin index for a given physical point p
virtual const unsigned int & getNx() const
Number of subchannels in the -x direction.
const ConsoleStream _console
virtual Real value(Real t, const Point &p) const
virtual const std::vector< unsigned int > & getChannelPins(unsigned int i_chan) const override
Return a vector of pin indices for a given channel index.
bool absoluteFuzzyGreaterThan(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
virtual const Real & getHeatedLengthEntry() const
Return unheated length at entry.
std::string _filename
The name of the radial power profile file.
virtual void initialSetup() override
registerMooseObjectRenamed("SubChannelApp", QuadPowerIC, "06/30/2025 24:00", SCMQuadPowerIC)
Sets the axial heat rate for each pin according to a radial power distribution and a user defined axi...