www.mooseframework.org
RichardsDensityVDW.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 #include "RichardsDensityVDW.h"
8 #include "libmesh/utility.h"
9 
10 template <>
11 InputParameters
13 {
14  InputParameters params = validParams<RichardsDensity>();
15  params.addRequiredRangeCheckedParam<Real>(
16  "a",
17  "a > 0",
18  "Parameter 'a' in the van der Waals expression (P + n^2 a/V^2)(V - nb) = nRT. "
19  "Example for methane 0.2303 Pa m^6 mol^-2");
20  params.addRequiredRangeCheckedParam<Real>(
21  "b",
22  "b > 0",
23  "Parameter 'b' in the van der Waals expression (P + n^2 a/V^2)(V - nb) = nRT. "
24  "Example for methane 4.31E-5 m^3/mol");
25  params.addRequiredRangeCheckedParam<Real>(
26  "temperature", "temperature > 0", "Temperature in Kelvin");
27  params.addRequiredRangeCheckedParam<Real>(
28  "molar_mass",
29  "molar_mass > 0",
30  "Molar mass of the gas. Example for methane 16.04246E-3 kg/mol");
31  params.addRangeCheckedParam<Real>("infinity_ratio",
32  10,
33  "infinity_ratio > 0",
34  "For P<0 the density is not physically defined, "
35  "but numerically it is advantageous to define "
36  "it: density(P=-infinity) = "
37  "-infinity_ratio*molar_mass, and density tends "
38  "exponentially towards this value as P -> "
39  "-infinity. (Units are mol/m^3).");
40  params.addClassDescription("Density of van der Waals gas.");
41  return params;
42 }
43 
44 RichardsDensityVDW::RichardsDensityVDW(const InputParameters & parameters)
45  : RichardsDensity(parameters),
46  _a(getParam<Real>("a")),
47  _b(getParam<Real>("b")),
48  _rt(getParam<Real>("temperature") * 8.314472), // multiply by gas constant
49  _molar_mass(getParam<Real>("molar_mass")),
50  _infinity_ratio(getParam<Real>("infinity_ratio")),
51  _rhs(_rt * _b / _a),
52  _b2oa(_b * _b / _a)
53 {
54  _vdw0 = densityVDW(0);
56 }
57 
58 Real
60 {
61  // transform (P + a/V^2)(V-b) = RT to
62  // (y + x^2)(1/x - 1) = rhs
63  // with: y = b^2*P/a, x = b/V, rhs = RT*b/a
64  // Then solve for x, using mathematica
65  // Then density = molar_mass/V = molar_mass*x/b
66  Real y = _b2oa * p;
67  Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
68  Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
69  Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
70  Real x = 1.0 / 3.0;
71  x += std::cbrt(2.0) * (-1.0 + 3.0 * _rhs + 3.0 * y) / (3.0 * cr);
72  x -= cr / (3.0 * std::cbrt(2.0));
73  return _molar_mass * x / _b;
74 }
75 
76 Real
78 {
79  if (p >= 0)
80  return densityVDW(p) - _vdw0;
81  else
82  return _infinity_ratio * _molar_mass * (std::exp(_slope0 * p) - 1.0);
83 }
84 
85 Real
87 {
88  if (p >= 0)
89  {
90  Real y = _b2oa * p;
91  Real dy = _b2oa;
92  Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
93  Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
94  Real dsq = 0.5 / sq * (4.0 * 3.0 * Utility::pow<2>(-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy +
95  2.0 * (-2.0 - 18.0 * y + 9.0 * _rhs) * (-18.0 * dy));
96  Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
97  Real dcr =
98  1.0 / 3.0 * std::pow(-2.0 + 9.0 * _rhs - 18.0 * y + sq, -2.0 / 3.0) * (-18.0 * dy + dsq);
99  Real dx = std::cbrt(2.0) *
100  ((3.0 * dy) / (3.0 * cr) + (-1.0 + 3.0 * _rhs + 3.0 * y) / 3.0 * (-dcr / (cr * cr)));
101  dx -= dcr / (3.0 * std::cbrt(2.0));
102  return _molar_mass * dx / _b;
103  }
104  else
105  return _infinity_ratio * _molar_mass * _slope0 * std::exp(_slope0 * p);
106 }
107 
108 Real
110 {
111  if (p >= 0)
112  {
113  Real y = _b2oa * p;
114  Real dy = _b2oa;
115  Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
116  Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
117  Real dsq = 0.5 / sq * (4.0 * 3.0 * Utility::pow<2>(-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy +
118  2.0 * (-2.0 - 18.0 * y + 9.0 * _rhs) * (-18.0 * dy));
119  Real d2sq = -dsq * dsq / sq;
120  d2sq += 0.5 / sq * (4.0 * 3.0 * 2.0 * (-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy * 3.0 * dy +
121  2.0 * (-18.0 * dy) * (-18.0 * dy));
122  Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
123  Real dcr =
124  1.0 / 3.0 * std::pow(-2.0 + 9.0 * _rhs - 18.0 * y + sq, -2.0 / 3.0) * (-18.0 * dy + dsq);
125  Real d2cr = 1.0 / 3.0 * (-2.0 / 3.0) * std::pow(-2 + 9 * _rhs - 18 * y + sq, -5. / 3.) *
126  Utility::pow<2>(-18 * dy + dsq);
127  d2cr += 1.0 / 3.0 * std::pow(-2 + 9 * _rhs - 18 * y + sq, -2. / 3.) * d2sq;
128  // Real dx = std::pow(2, 1.0/3.0)*( (3*dy)/3/cr + (-1 + 3*_rhs + 3*y)/3*(-dcr/cr/cr));
129  // dx -= dcr/3/std::pow(2, 1.0/3.0);
130  Real d2x = std::cbrt(2.0) *
131  (-(3.0 * dy) * dcr / (3.0 * cr * cr) + 3.0 * dy / 3.0 * (-dcr / (cr * cr)) +
132  (-1.0 + 3.0 * _rhs + 3.0 * y) / 3.0 *
133  (-d2cr / (cr * cr) + 2.0 * dcr * dcr / Utility::pow<3>(cr)));
134  d2x -= d2cr / (3.0 * std::cbrt(2.0));
135  return _molar_mass * d2x / _b;
136  }
137  else
138  return _infinity_ratio * _molar_mass * Utility::pow<2>(_slope0) * std::exp(_slope0 * p);
139 }
RichardsDensityVDW(const InputParameters &parameters)
Real _infinity_ratio
density at P=-infinity is _infinity_ratio*_molar_mass
Real _b
van der Waals b
Real _slope0
(1/_molar_mass)*d(density)/dP at P=0
Real d2density(Real p) const
second derivative of fluid density wrt porepressure
Real densityVDW(Real p) const
Density according to the van der Waals expression This is modified to yield density(p) as noted above...
Real _vdw0
density at P=0 according to the van der Waals expression
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
Real ddensity(Real p) const
derivative of fluid density wrt porepressure
Real density(Real p) const
fluid density as a function of porepressure
InputParameters validParams< RichardsDensityVDW >()
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
Real _molar_mass
molar mass of gas
InputParameters validParams< RichardsDensity >()