www.mooseframework.org
TabulatedFluidProperties.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
9 #include "BicubicSplineInterpolation.h"
10 #include "MooseUtils.h"
11 
12 // C++ includes
13 #include <fstream>
14 #include <ctime>
15 
16 template <>
17 InputParameters
19 {
20  InputParameters params = validParams<SinglePhaseFluidPropertiesPT>();
21  params.addParam<FileName>("fluid_property_file",
22  "fluid_properties.csv",
23  "Name of the csv file containing the tabulated fluid property data");
24  params.addRangeCheckedParam<Real>("temperature_min",
25  300.0,
26  "temperature_min > 0",
27  "Minimum temperature for tabulated data. Default is 300 K)");
28  params.addParam<Real>(
29  "temperature_max", 500.0, "Maximum temperature for tabulated data. Default is 500 K");
30  params.addRangeCheckedParam<Real>("pressure_min",
31  1.0e5,
32  "pressure_min > 0",
33  "Minimum pressure for tabulated data. Default is 0.1 MPa)");
34  params.addParam<Real>(
35  "pressure_max", 50.0e6, "Maximum pressure for tabulated data. Default is 50 MPa");
36  params.addRangeCheckedParam<unsigned int>(
37  "num_T", 100, "num_T > 0", "Number of points to divide temperature range. Default is 100");
38  params.addRangeCheckedParam<unsigned int>(
39  "num_p", 100, "num_p > 0", "Number of points to divide pressure range. Default is 100");
40  params.addRequiredParam<UserObjectName>("fp", "The name of the FluidProperties UserObject");
41  params.addClassDescription(
42  "Fluid properties using bicubic spline interpolation on tabulated values provided");
43  return params;
44 }
45 
46 TabulatedFluidProperties::TabulatedFluidProperties(const InputParameters & parameters)
47  : SinglePhaseFluidPropertiesPT(parameters),
48  _file_name(getParam<FileName>("fluid_property_file")),
49  _temperature_min(getParam<Real>("temperature_min")),
50  _temperature_max(getParam<Real>("temperature_max")),
51  _pressure_min(getParam<Real>("pressure_min")),
52  _pressure_max(getParam<Real>("pressure_max")),
53  _num_T(getParam<unsigned int>("num_T")),
54  _num_p(getParam<unsigned int>("num_p")),
55  _fp(getUserObject<SinglePhaseFluidPropertiesPT>("fp")),
56  _csv_reader(_file_name, &_communicator)
57 {
58  // Sanity check on minimum and maximum temperatures and pressures
60  mooseError("temperature_max must be greater than temperature_min in ", name());
62  mooseError("pressure_max must be greater than pressure_min in ", name());
63 
64  // Lines starting with # are treated as comments
65  _csv_reader.setComment("#");
66 }
67 
69 
70 void
72 {
73  // Check to see if _file_name supplied exists. If it does, that data
74  // will be used. If it does not exist, data will be generated and then
75  // written to _file_name.
76  std::ifstream file(_file_name.c_str());
77  if (file.good())
78  {
79  _console << "Reading tabulated properties from " << _file_name << "\n";
80  _csv_reader.read();
81 
82  const std::vector<std::string> & column_names = _csv_reader.getNames();
83 
84  // Check that all required columns are present
85  for (std::size_t i = 0; i < _required_columns.size(); ++i)
86  {
87  if (std::find(column_names.begin(), column_names.end(), _required_columns[i]) ==
88  column_names.end())
89  mooseError("No ",
91  " data read in ",
92  _file_name,
93  ". A column named ",
95  " must be present");
96  }
97 
98  std::map<std::string, unsigned int> data_index;
99  for (std::size_t i = 0; i < _required_columns.size(); ++i)
100  {
101  auto it = std::find(column_names.begin(), column_names.end(), _required_columns[i]);
102  data_index[_required_columns[i]] = std::distance(column_names.begin(), it);
103  }
104 
105  const std::vector<std::vector<Real>> & column_data = _csv_reader.getData();
106 
107  // Extract the pressure and temperature data vectors
108  _pressure = column_data[data_index.find("pressure")->second];
109  _temperature = column_data[data_index.find("temperature")->second];
110 
111  // Pressure and temperature data contains duplicates due to the csv format.
112  // First, check that pressure is monotonically increasing
113  if (!std::is_sorted(_pressure.begin(), _pressure.end()))
114  mooseError("The column data for pressure is not monotonically increasing in ", _file_name);
115 
116  // The first pressure value is repeated for each temperature value. Counting the
117  // number of repeats provides the number of temperature values
118  auto num_T = std::count(_pressure.begin(), _pressure.end(), _pressure.front());
119 
120  // Now remove the duplicates in the pressure vector
121  auto last_unique = std::unique(_pressure.begin(), _pressure.end());
122  _pressure.erase(last_unique, _pressure.end());
123  _num_p = _pressure.size();
124 
125  // Check that the number of rows in the csv file is equal to _num_p * _num_T
126  if (column_data[0].size() != _num_p * static_cast<unsigned int>(num_T))
127  mooseError("The number of rows in ",
128  _file_name,
129  " is not equal to the number of unique pressure values ",
130  _num_p,
131  " multiplied by the number of unique temperature values ",
132  num_T);
133 
134  // Need to make sure that the temperature values are provided in ascending order
135  // as well as duplicated for each pressure value
136  std::vector<Real> temp0(_temperature.begin(), _temperature.begin() + num_T);
137  if (!std::is_sorted(temp0.begin(), temp0.end()))
138  mooseError("The column data for temperature is not monotonically increasing in ", _file_name);
139 
140  auto it_temp = _temperature.begin() + num_T;
141  for (std::size_t i = 1; i < _pressure.size(); ++i)
142  {
143  std::vector<Real> temp(it_temp, it_temp + num_T);
144  if (temp != temp0)
145  mooseError("Temperature values for pressure ",
146  _pressure[i],
147  " are not identical to values for ",
148  _pressure[0]);
149 
150  std::advance(it_temp, num_T);
151  }
152 
153  // At this point, all temperature data has been provided in ascending order
154  // identically for each pressure value, so we can just keep the first range
155  _temperature.erase(_temperature.begin() + num_T, _temperature.end());
156  _num_T = _temperature.size();
157 
158  // Minimum and maximum pressure and temperature. Note that _pressure and
159  // _temperature are sorted
160  _pressure_min = _pressure.front();
161  _pressure_max = _pressure.back();
162  _temperature_min = _temperature.front();
164 
165  // Extract the fluid property data and reshape into 2D arrays for interpolation
166  reshapeData2D(_num_p, _num_T, column_data[data_index.find("density")->second], _density);
167  reshapeData2D(_num_p, _num_T, column_data[data_index.find("enthalpy")->second], _enthalpy);
169  _num_p, _num_T, column_data[data_index.find("internal_energy")->second], _internal_energy);
170  }
171  else
172  {
173  _console << "No tabulated properties file named " << _file_name << " exists.\n"
174  << "Generating tabulated data and writing output to " << _file_name << "\n";
175 
177 
178  // Write tabulated data to file
179  writeTabulatedData(_file_name);
180  }
181 
182  // Construct bicubic splines from tabulated data
183  _density_ipol = libmesh_make_unique<BicubicSplineInterpolation>();
184  _internal_energy_ipol = libmesh_make_unique<BicubicSplineInterpolation>();
185  _enthalpy_ipol = libmesh_make_unique<BicubicSplineInterpolation>();
186 
187  _density_ipol->setData(
189 
190  _internal_energy_ipol->setData(
192 
193  _enthalpy_ipol->setData(
195 }
196 
197 std::string
199 {
200  return _fp.fluidName();
201 }
202 
203 Real
205 {
206  return _fp.molarMass();
207 }
208 
209 Real
211 {
212  checkInputVariables(pressure, temperature);
213  return _density_ipol->sample(pressure, temperature);
214 }
215 
216 void
218  Real pressure, Real temperature, Real & rho, Real & drho_dp, Real & drho_dT) const
219 {
220  checkInputVariables(pressure, temperature);
221  rho = _density_ipol->sample(pressure, temperature);
222  drho_dp = _density_ipol->sampleDerivative(pressure, temperature, _wrt_p);
223  drho_dT = _density_ipol->sampleDerivative(pressure, temperature, _wrt_T);
224 }
225 
226 Real
228 {
229  checkInputVariables(pressure, temperature);
230  return _internal_energy_ipol->sample(pressure, temperature);
231 }
232 
233 void
235  Real pressure, Real temperature, Real & e, Real & de_dp, Real & de_dT) const
236 {
237  checkInputVariables(pressure, temperature);
238  e = _internal_energy_ipol->sample(pressure, temperature);
239  de_dp = _internal_energy_ipol->sampleDerivative(pressure, temperature, _wrt_p);
240  de_dT = _internal_energy_ipol->sampleDerivative(pressure, temperature, _wrt_T);
241 }
242 
243 void
245  Real temperature,
246  Real & rho,
247  Real & drho_dp,
248  Real & drho_dT,
249  Real & e,
250  Real & de_dp,
251  Real & de_dT) const
252 {
253  checkInputVariables(pressure, temperature);
254  rho_dpT(pressure, temperature, rho, drho_dp, drho_dT);
255  e_dpT(pressure, temperature, e, de_dp, de_dT);
256 }
257 
258 Real
260 {
261  checkInputVariables(pressure, temperature);
262  return _enthalpy_ipol->sample(pressure, temperature);
263 }
264 
265 void
267  Real pressure, Real temperature, Real & h, Real & dh_dp, Real & dh_dT) const
268 {
269  checkInputVariables(pressure, temperature);
270  h = _enthalpy_ipol->sample(pressure, temperature);
271  dh_dp = _enthalpy_ipol->sampleDerivative(pressure, temperature, _wrt_p);
272  dh_dT = _enthalpy_ipol->sampleDerivative(pressure, temperature, _wrt_T);
273 }
274 
275 Real
277 {
278  Real rho = this->rho(pressure, temperature);
279  return this->mu_from_rho_T(rho, temperature);
280 }
281 
282 void
284  Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const
285 {
286  Real rho, drho_dp, drho_dT;
287  this->rho_dpT(pressure, temperature, rho, drho_dp, drho_dT);
288  Real dmu_drho;
289  this->mu_drhoT_from_rho_T(rho, temperature, drho_dT, mu, dmu_drho, dmu_dT);
290  dmu_dp = dmu_drho * drho_dp;
291 }
292 
293 Real
295 {
296  return _fp.mu_from_rho_T(density, temperature);
297 }
298 
299 void
301  Real temperature,
302  Real ddensity_dT,
303  Real & mu,
304  Real & dmu_drho,
305  Real & dmu_dT) const
306 {
307  _fp.mu_drhoT_from_rho_T(density, temperature, ddensity_dT, mu, dmu_drho, dmu_dT);
308 }
309 
310 Real
312 {
313  return _fp.c(pressure, temperature);
314 }
315 
316 Real
318 {
319  return _fp.cp(pressure, temperature);
320 }
321 
322 Real
324 {
325  return _fp.cv(pressure, temperature);
326 }
327 
328 Real
330 {
331  Real rho = this->rho(pressure, temperature);
332  return this->k_from_rho_T(rho, temperature);
333 }
334 
335 void
337  Real /*pressure*/, Real /*temperature*/, Real & /*k*/, Real & /*dk_dp*/, Real & /*dk_dT*/) const
338 {
339  mooseError(name(), "k_dpT() is not implemented");
340 }
341 
342 Real
344 {
345  return _fp.k_from_rho_T(density, temperature);
346 }
347 
348 Real
350 {
351  return _fp.s(pressure, temperature);
352 }
353 
354 Real
356 {
357  return _fp.beta(pressure, temperature);
358 }
359 
360 Real
362 {
363  return _fp.henryConstant(temperature);
364 }
365 
366 void
367 TabulatedFluidProperties::henryConstant_dT(Real temperature, Real & Kh, Real & dKh_dT) const
368 {
369  _fp.henryConstant_dT(temperature, Kh, dKh_dT);
370 }
371 
372 void
374 {
375  if (processor_id() == 0)
376  {
377  MooseUtils::checkFileWriteable(file_name);
378 
379  std::ofstream file_out(file_name.c_str());
380 
381  // Write out date and fluid type
382  time_t now = time(&now);
383  file_out << "# " << _fp.fluidName() << " properties created by TabulatedFluidProperties on "
384  << ctime(&now) << "\n";
385 
386  std::string column_names{"pressure, temperature, density, enthalpy, internal_energy"};
387 
388  // Write out column names
389  file_out << column_names << "\n";
390 
391  // Write out the fluid property data
392  for (unsigned int i = 0; i < _num_p; ++i)
393  for (unsigned int j = 0; j < _num_T; ++j)
394  file_out << _pressure[i] << ", " << _temperature[j] << ", " << _density[i][j] << ", "
395  << _enthalpy[i][j] << ", " << _internal_energy[i][j] << "\n";
396  }
397 }
398 
399 void
401 {
402  _pressure.resize(_num_p);
403  _temperature.resize(_num_T);
404 
405  _density.resize(_num_p);
406  _internal_energy.resize(_num_p);
407  _enthalpy.resize(_num_p);
408 
409  for (unsigned int i = 0; i < _num_p; ++i)
410  {
411  _density[i].resize(_num_T);
412  _internal_energy[i].resize(_num_T);
413  _enthalpy[i].resize(_num_T);
414  }
415 
416  // Temperature is divided equally into _num_T segments
417  Real delta_T = (_temperature_max - _temperature_min) / static_cast<Real>(_num_T - 1);
418 
419  for (unsigned int j = 0; j < _num_T; ++j)
420  _temperature[j] = _temperature_min + j * delta_T;
421 
422  // Divide the pressure into _num_p equal segments
423  Real delta_p = (_pressure_max - _pressure_min) / static_cast<Real>(_num_p - 1);
424 
425  for (unsigned int i = 0; i < _num_p; ++i)
426  _pressure[i] = _pressure_min + i * delta_p;
427 
428  // Generate the tabulated data at the pressure and temperature points
429  for (unsigned int i = 0; i < _num_p; ++i)
430  for (unsigned int j = 0; j < _num_T; ++j)
431  {
432  _density[i][j] = _fp.rho(_pressure[i], _temperature[j]);
433  _internal_energy[i][j] = _fp.e(_pressure[i], _temperature[j]);
434  _enthalpy[i][j] = _fp.h(_pressure[i], _temperature[j]);
435  }
436 }
437 
438 void
440  unsigned int ncol,
441  const std::vector<Real> & vec,
442  std::vector<std::vector<Real>> & mat)
443 {
444  if (!vec.empty())
445  {
446  mat.resize(nrow);
447 
448  for (unsigned int i = 0; i < nrow; ++i)
449  for (unsigned int j = 0; j < ncol; ++j)
450  mat[i].push_back(vec[i * ncol + j]);
451  }
452 }
453 
454 void
456 {
457  if (pressure < _pressure_min || pressure > _pressure_max)
458  mooseError("Pressure ",
459  pressure,
460  " is outside the range of tabulated pressure (",
462  ", ",
463  _pressure_max,
464  ".");
465 
466  if (temperature < _temperature_min || temperature > _temperature_max)
467  mooseError("Temperature ",
468  temperature,
469  " is outside the range of tabulated temperature (",
471  ", ",
472  _temperature_max,
473  ".");
474 }
virtual void initialSetup() override
std::unique_ptr< BicubicSplineInterpolation > _enthalpy_ipol
Interpoled enthalpy.
virtual Real beta(Real pressure, Real temperature) const =0
Thermal expansion coefficient.
virtual std::string fluidName() const =0
Fluid name.
virtual Real k(Real pressure, Real temperature) const override
Thermal conductivity.
virtual Real c(Real pressure, Real temperature) const override
Speed of sound.
virtual Real henryConstant(Real temperature) const override
Henry&#39;s law constant for dissolution in water.
std::vector< std::vector< Real > > _internal_energy
Tabulated internal energy.
virtual Real cp(Real pressure, Real temperature) const =0
Isobaric specific heat capacity.
Real _temperature_max
Maximum temperature in tabulated data.
const std::vector< std::string > _required_columns
List of required column names to be read.
virtual Real rho(Real pressure, Real temperature) const override
Density.
std::unique_ptr< BicubicSplineInterpolation > _internal_energy_ipol
Interpoled internal energy.
const unsigned int _wrt_T
Index for derivatives wrt temperature.
virtual Real mu_from_rho_T(Real density, Real temperature) const override
const std::string density
Definition: NS.h:15
Real _pressure_max
Maximum pressure in tabulated data.
std::vector< std::vector< Real > > _density
Tabulated density.
virtual void mu_dpT(Real pressure, Real temperature, Real &mu, Real &dmu_dp, Real &dmu_dT) const override
MooseUtils::DelimitedFileReader _csv_reader
The MOOSE delimited file reader.
virtual void henryConstant_dT(Real temperature, Real &Kh, Real &dKh_dT) const override
Henry&#39;s law constant for dissolution in water and derivative wrt temperature.
virtual Real c(Real pressure, Real temperature) const =0
Speed of sound.
virtual Real mu_from_rho_T(Real density, Real temperature) const =0
const std::string temperature
Definition: NS.h:25
virtual Real beta(Real pressure, Real temperature) const override
Thermal expansion coefficient.
void checkInputVariables(Real pressure, Real temperature) const
Checks that the inputs are within the range of the tabulated data, and throws an error if they are no...
virtual std::string fluidName() const override
Fluid name.
virtual Real s(Real pressure, Real temperature) const override
Specific entropy.
Common class for single phase fluid properties using a pressure and temperature formulation.
const unsigned int _wrt_p
Index for derivatives wrt pressure.
virtual Real henryConstant(Real temperature) const =0
Henry&#39;s law constant for dissolution in water.
virtual Real molarMass() const =0
Molar mass.
InputParameters validParams< SinglePhaseFluidPropertiesPT >()
virtual void henryConstant_dT(Real temperature, Real &Kh, Real &dKh_dT) const =0
Henry&#39;s law constant for dissolution in water and derivative wrt temperature.
Real _temperature_min
Minimum temperature in tabulated data.
std::unique_ptr< BicubicSplineInterpolation > _density_ipol
Interpoled density.
virtual void e_dpT(Real pressure, Real temperature, Real &e, Real &de_dp, Real &de_dT) const override
Internal energy and its derivatives wrt pressure and temperature.
virtual void h_dpT(Real pressure, Real temperature, Real &h, Real &dh_dp, Real &dh_dT) const override
Enthalpy and its derivatives wrt pressure and temperature.
virtual void mu_drhoT_from_rho_T(Real density, Real temperature, Real ddensity_dT, Real &mu, Real &dmu_drho, Real &dmu_dT) const =0
Dynamic viscosity and its derivatives wrt density and temperature.
const SinglePhaseFluidPropertiesPT & _fp
SinglePhaseFluidPropertiesPT UserObject.
std::vector< std::vector< Real > > _enthalpy
Tabulated enthalpy.
virtual void rho_e_dpT(Real pressure, Real temperature, Real &rho, Real &drho_dp, Real &drho_dT, Real &e, Real &de_dp, Real &de_dT) const override
Density and internal energy and their derivatives wrt pressure and temperature.
virtual Real e(Real pressure, Real temperature) const =0
Internal energy.
TabulatedFluidProperties(const InputParameters &parameters)
virtual void rho_dpT(Real pressure, Real temperature, Real &rho, Real &drho_dp, Real &drho_dT) const override
Density and its derivatives wrt pressure and temperature.
virtual Real mu(Real pressure, Real temperature) const override
std::vector< Real > _pressure
Pressure vector.
std::vector< Real > _temperature
Temperature vector.
virtual Real cv(Real pressure, Real temperature) const override
Isochoric specific heat.
virtual void k_dpT(Real pressure, Real temperature, Real &k, Real &dk_dp, Real &dk_dT) const override
Thermal conductivity and its derivatives wrt pressure and temperature.
unsigned int _num_T
Number of temperature points in the tabulated data.
Real _pressure_min
Minimum pressure in tabulated data.
unsigned int _num_p
Number of pressure points in the tabulated data.
virtual Real h(Real p, Real T) const =0
Specific enthalpy.
virtual Real cv(Real pressure, Real temperature) const =0
Isochoric specific heat.
virtual Real rho(Real pressure, Real temperature) const =0
Density.
const std::string pressure
Definition: NS.h:24
virtual Real e(Real pressure, Real temperature) const override
Internal energy.
void reshapeData2D(unsigned int nrow, unsigned int ncol, const std::vector< Real > &vec, std::vector< std::vector< Real >> &mat)
Forms a 2D matrix from a single std::vector.
virtual Real k_from_rho_T(Real density, Real temperature) const =0
Thermal conductivity as a function of density and temperature.
virtual void mu_drhoT_from_rho_T(Real density, Real temperature, Real ddensity_dT, Real &mu, Real &dmu_drho, Real &dmu_dT) const override
Dynamic viscosity and its derivatives wrt density and temperature.
virtual Real h(Real p, Real T) const override
Specific enthalpy.
std::vector< Real > _drho_dp_0
Derivatives along the boundary.
virtual Real s(Real pressure, Real temperature) const =0
Specific entropy.
virtual Real cp(Real pressure, Real temperature) const override
Isobaric specific heat capacity.
virtual void generateTabulatedData()
Generates a table of fluid properties by looping over pressure and temperature and calculating proper...
void writeTabulatedData(std::string file_name)
Writes tabulated data to a file.
virtual Real k_from_rho_T(Real density, Real temperature) const override
Thermal conductivity as a function of density and temperature.
virtual Real molarMass() const override
Molar mass.
InputParameters validParams< TabulatedFluidProperties >()
FileName _file_name
File name of tabulated data file.