www.mooseframework.org
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
RichardsDensityVDW Class Reference

Density of a gas according to the van der Waals expression (P + n^2 a/V^2)(V - nb) = nRT How density is calculated: given P, (1/V) is calculated for n=1 and rho = molar_mass*(1/V). More...

#include <RichardsDensityVDW.h>

Inheritance diagram for RichardsDensityVDW:
[legend]

Public Member Functions

 RichardsDensityVDW (const InputParameters &parameters)
 
Real density (Real p) const
 fluid density as a function of porepressure More...
 
Real ddensity (Real p) const
 derivative of fluid density wrt porepressure More...
 
Real d2density (Real p) const
 second derivative of fluid density wrt porepressure More...
 
void initialize ()
 
void execute ()
 
void finalize ()
 

Protected Member Functions

Real densityVDW (Real p) const
 Density according to the van der Waals expression This is modified to yield density(p) as noted above. More...
 

Protected Attributes

Real _a
 van der Waals a More...
 
Real _b
 van der Waals b More...
 
Real _rt
 R*T (gas constant * temperature) More...
 
Real _molar_mass
 molar mass of gas More...
 
Real _infinity_ratio
 density at P=-infinity is _infinity_ratio*_molar_mass More...
 
Real _rhs
 R*T*b/a. More...
 
Real _b2oa
 b^2/a More...
 
Real _vdw0
 density at P=0 according to the van der Waals expression More...
 
Real _slope0
 (1/_molar_mass)*d(density)/dP at P=0 More...
 

Detailed Description

Density of a gas according to the van der Waals expression (P + n^2 a/V^2)(V - nb) = nRT How density is calculated: given P, (1/V) is calculated for n=1 and rho = molar_mass*(1/V).

The density is actually modified from this rho in the following ways: (1) density = rho - rho(P=0). This is so that density(P=0)=0, which is physically correct, but the wan der Waals expression does not hold close to a vacuum. However, it can be vital that density(P=0)=0 numerically, otherwise fluid can be withdrawn from a node where there shouldn't be any fluid at P=0. This is a miniscule correction, of many orders of magnitude below the precision of a, b, R or T (2) For P<0 we enter an unphysical region, but this may be sampled by the Newton process as the algorithm iterates towards the solution. Therefore i set density = infinity_ratio*molar_mass*(exp(-c*P) - 1) in this region where c is chosen so the derivatives match at P=0

Definition at line 35 of file RichardsDensityVDW.h.

Constructor & Destructor Documentation

RichardsDensityVDW::RichardsDensityVDW ( const InputParameters &  parameters)

Definition at line 44 of file RichardsDensityVDW.C.

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 }
RichardsDensity(const InputParameters &parameters)
Real _infinity_ratio
density at P=-infinity is _infinity_ratio*_molar_mass
Real _b
van der Waals b
Real _a
van der Waals a
Real _slope0
(1/_molar_mass)*d(density)/dP at P=0
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
Real _rt
R*T (gas constant * temperature)
Real ddensity(Real p) const
derivative of fluid density wrt porepressure
Real _molar_mass
molar mass of gas

Member Function Documentation

Real RichardsDensityVDW::d2density ( Real  p) const
virtual

second derivative of fluid density wrt porepressure

Parameters
pporepressure

Implements RichardsDensity.

Definition at line 109 of file RichardsDensityVDW.C.

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 }
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
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
Real _molar_mass
molar mass of gas
Real RichardsDensityVDW::ddensity ( Real  p) const
virtual

derivative of fluid density wrt porepressure

Parameters
pporepressure

Implements RichardsDensity.

Definition at line 86 of file RichardsDensityVDW.C.

Referenced by RichardsDensityVDW().

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 }
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
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
Real _molar_mass
molar mass of gas
Real RichardsDensityVDW::density ( Real  p) const
virtual

fluid density as a function of porepressure

Parameters
pporepressure

Implements RichardsDensity.

Definition at line 77 of file RichardsDensityVDW.C.

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 }
Real _infinity_ratio
density at P=-infinity is _infinity_ratio*_molar_mass
Real _slope0
(1/_molar_mass)*d(density)/dP at P=0
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
Real _molar_mass
molar mass of gas
Real RichardsDensityVDW::densityVDW ( Real  p) const
protected

Density according to the van der Waals expression This is modified to yield density(p) as noted above.

Parameters
pporepressure

Definition at line 59 of file RichardsDensityVDW.C.

Referenced by density(), and RichardsDensityVDW().

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 }
Real _b
van der Waals b
Real _molar_mass
molar mass of gas
void RichardsDensity::execute ( )
inherited

Definition at line 32 of file RichardsDensity.C.

33 {
34 }
void RichardsDensity::finalize ( )
inherited

Definition at line 37 of file RichardsDensity.C.

38 {
39 }
void RichardsDensity::initialize ( )
inherited

Definition at line 27 of file RichardsDensity.C.

28 {
29 }

Member Data Documentation

Real RichardsDensityVDW::_a
protected

van der Waals a

Definition at line 60 of file RichardsDensityVDW.h.

Real RichardsDensityVDW::_b
protected

van der Waals b

Definition at line 63 of file RichardsDensityVDW.h.

Referenced by d2density(), ddensity(), and densityVDW().

Real RichardsDensityVDW::_b2oa
protected

b^2/a

Definition at line 78 of file RichardsDensityVDW.h.

Referenced by d2density(), ddensity(), and densityVDW().

Real RichardsDensityVDW::_infinity_ratio
protected

density at P=-infinity is _infinity_ratio*_molar_mass

Definition at line 72 of file RichardsDensityVDW.h.

Referenced by d2density(), ddensity(), density(), and RichardsDensityVDW().

Real RichardsDensityVDW::_molar_mass
protected

molar mass of gas

Definition at line 69 of file RichardsDensityVDW.h.

Referenced by d2density(), ddensity(), density(), densityVDW(), and RichardsDensityVDW().

Real RichardsDensityVDW::_rhs
protected

R*T*b/a.

Definition at line 75 of file RichardsDensityVDW.h.

Referenced by d2density(), ddensity(), and densityVDW().

Real RichardsDensityVDW::_rt
protected

R*T (gas constant * temperature)

Definition at line 66 of file RichardsDensityVDW.h.

Real RichardsDensityVDW::_slope0
protected

(1/_molar_mass)*d(density)/dP at P=0

Definition at line 84 of file RichardsDensityVDW.h.

Referenced by d2density(), ddensity(), density(), and RichardsDensityVDW().

Real RichardsDensityVDW::_vdw0
protected

density at P=0 according to the van der Waals expression

Definition at line 81 of file RichardsDensityVDW.h.

Referenced by density(), and RichardsDensityVDW().


The documentation for this class was generated from the following files: