www.mooseframework.org
RichardsDensityConstBulkCut.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 // Fluid density assuming constant bulk modulus, but reducing via a cubic to zero at p=zero_point
11 //
13 
15 
18 {
20  params.addRequiredParam<Real>("dens0", "Reference density of fluid. Eg 1000");
21  params.addRequiredParam<Real>("bulk_mod", "Bulk modulus of fluid. Eg 2E9");
22  params.addParam<Real>("cut_limit",
23  1E5,
24  "For porepressure > cut_limit, density = dens0*Exp(pressure/bulk). For "
25  "porepressure < cut_limie, density = cubic*dens0*Exp(pressure/bulk), where "
26  "cubic=1 for pressure=cut_limit, and cubic=0 for pressure<=zero_point");
27  params.addParam<Real>("zero_point",
28  0,
29  "For porepressure > cut_limit, density = dens0*Exp(pressure/bulk). For "
30  "porepressure < cut_limie, density = cubic*dens0*Exp(pressure/bulk), where "
31  "cubic=1 for pressure=cut_limit, and cubic=0 for pressure<=zero_point");
32  params.addClassDescription(
33  "Fluid density assuming constant bulk modulus. dens0 * Exp(pressure/bulk)");
34  return params;
35 }
36 
38  : RichardsDensity(parameters),
39  _dens0(getParam<Real>("dens0")),
40  _bulk(getParam<Real>("bulk_mod")),
41  _cut_limit(getParam<Real>("cut_limit")),
42  _zero_point(getParam<Real>("zero_point")),
43  _c3(std::pow(_cut_limit - _zero_point, 3))
44 {
45  if (_zero_point >= _cut_limit)
46  mooseError("RichardsDensityConstantbulkCut: zero_point must be less than cut_limit");
47 }
48 
49 Real
51 {
52  if (p <= _zero_point)
53  return 0;
54 
55  Real unmodified = _dens0 * std::exp(p / _bulk);
56  if (p >= _cut_limit)
57  return unmodified;
58 
59  Real modifier =
60  (3 * _cut_limit - 2 * p - _zero_point) * (p - _zero_point) * (p - _zero_point) / _c3;
61  return modifier * unmodified;
62 }
63 
64 Real
66 {
67  if (p <= _zero_point)
68  return 0;
69 
70  Real unmodified = _dens0 * std::exp(p / _bulk);
71  if (p >= _cut_limit)
72  return unmodified / _bulk;
73 
74  Real modifier =
75  (3 * _cut_limit - 2 * p - _zero_point) * (p - _zero_point) * (p - _zero_point) / _c3;
76  Real dmodifier = 6 * (_cut_limit - p) * (p - _zero_point) / _c3;
77  return (modifier / _bulk + dmodifier) * unmodified;
78 }
79 
80 Real
82 {
83  if (p <= _zero_point)
84  return 0;
85 
86  Real unmodified = _dens0 * std::exp(p / _bulk);
87  if (p >= _cut_limit)
88  return unmodified / _bulk / _bulk;
89 
90  Real modifier =
91  (3 * _cut_limit - 2 * p - _zero_point) * (p - _zero_point) * (p - _zero_point) / _c3;
92  Real dmodifier = 6 * (_cut_limit - p) * (p - _zero_point) / _c3;
93  Real d2modifier = (6 * _cut_limit + 6 * _zero_point - 12 * p) / _c3;
94  return (modifier / _bulk / _bulk + 2 * dmodifier / _bulk + d2modifier) * unmodified;
95 }
registerMooseObject("RichardsApp", RichardsDensityConstBulkCut)
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
Real _cut_limit
where the cubic starts
Real d2density(Real p) const
second derivative of fluid density wrt porepressure
static InputParameters validParams()
static InputParameters validParams()
void addRequiredParam(const std::string &name, const std::string &doc_string)
Real _dens0
density = _dens0*exp(p/_bulk), modified by cubic
Real _bulk
density = _dens0*exp(p/_bulk), modified by cubic
Real density(Real p) const
fluid density as a function of porepressure
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real _c3
(cut_limit-zero_point)^3
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
Fluid density assuming constant bulk modulus, for p>cut_limit Then following a cubic for zero_point <...
Real ddensity(Real p) const
derivative of fluid density wrt porepressure
Real _zero_point
where the density is zero
RichardsDensityConstBulkCut(const InputParameters &parameters)