RichardsDensityConstBulkCut.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
8 // Fluid density assuming constant bulk modulus, but reducing via a cubic to zero at p=zero_point
9 //
11
12 template <>
13 InputParameters
15 {
16  InputParameters params = validParams<RichardsDensity>();
17  params.addRequiredParam<Real>("dens0", "Reference density of fluid. Eg 1000");
18  params.addRequiredParam<Real>("bulk_mod", "Bulk modulus of fluid. Eg 2E9");
20  1E5,
21  "For porepressure > cut_limit, density = dens0*Exp(pressure/bulk). For "
22  "porepressure < cut_limie, density = cubic*dens0*Exp(pressure/bulk), where "
23  "cubic=1 for pressure=cut_limit, and cubic=0 for pressure<=zero_point");
25  0,
26  "For porepressure > cut_limit, density = dens0*Exp(pressure/bulk). For "
27  "porepressure < cut_limie, density = cubic*dens0*Exp(pressure/bulk), where "
28  "cubic=1 for pressure=cut_limit, and cubic=0 for pressure<=zero_point");
30  "Fluid density assuming constant bulk modulus. dens0 * Exp(pressure/bulk)");
31  return params;
32 }
33
35  : RichardsDensity(parameters),
36  _dens0(getParam<Real>("dens0")),
37  _bulk(getParam<Real>("bulk_mod")),
38  _cut_limit(getParam<Real>("cut_limit")),
39  _zero_point(getParam<Real>("zero_point")),
40  _c3(std::pow(_cut_limit - _zero_point, 3))
41 {
42  if (_zero_point >= _cut_limit)
43  mooseError("RichardsDensityConstantbulkCut: zero_point must be less than cut_limit");
44 }
45
46 Real
48 {
49  if (p <= _zero_point)
50  return 0;
51
52  Real unmodified = _dens0 * std::exp(p / _bulk);
53  if (p >= _cut_limit)
54  return unmodified;
55
56  Real modifier =
57  (3 * _cut_limit - 2 * p - _zero_point) * (p - _zero_point) * (p - _zero_point) / _c3;
58  return modifier * unmodified;
59 }
60
61 Real
63 {
64  if (p <= _zero_point)
65  return 0;
66
67  Real unmodified = _dens0 * std::exp(p / _bulk);
68  if (p >= _cut_limit)
69  return unmodified / _bulk;
70
71  Real modifier =
72  (3 * _cut_limit - 2 * p - _zero_point) * (p - _zero_point) * (p - _zero_point) / _c3;
73  Real dmodifier = 6 * (_cut_limit - p) * (p - _zero_point) / _c3;
74  return (modifier / _bulk + dmodifier) * unmodified;
75 }
76
77 Real
79 {
80  if (p <= _zero_point)
81  return 0;
82
83  Real unmodified = _dens0 * std::exp(p / _bulk);
84  if (p >= _cut_limit)
85  return unmodified / _bulk / _bulk;
86
87  Real modifier =
88  (3 * _cut_limit - 2 * p - _zero_point) * (p - _zero_point) * (p - _zero_point) / _c3;
89  Real dmodifier = 6 * (_cut_limit - p) * (p - _zero_point) / _c3;
90  Real d2modifier = (6 * _cut_limit + 6 * _zero_point - 12 * p) / _c3;
91  return (modifier / _bulk / _bulk + 2 * dmodifier / _bulk + d2modifier) * unmodified;
92 }
Real density(Real p) const
fluid density as a function of porepressure
Real _cut_limit
where the cubic starts
Real ddensity(Real p) const
derivative of fluid density wrt porepressure
Real _dens0
density = _dens0*exp(p/_bulk), modified by cubic
Real _bulk
density = _dens0*exp(p/_bulk), modified by cubic
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
InputParameters validParams< RichardsDensityConstBulkCut >()
Real _c3
(cut_limit-zero_point)^3
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
Real d2density(Real p) const
second derivative of fluid density wrt porepressure
InputParameters validParams< RichardsDensity >()
Real _zero_point
where the density is zero
RichardsDensityConstBulkCut(const InputParameters &parameters)