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 37 of file RichardsDensityVDW.h.

Constructor & Destructor Documentation

◆ RichardsDensityVDW()

RichardsDensityVDW::RichardsDensityVDW ( const InputParameters &  parameters)

Definition at line 49 of file RichardsDensityVDW.C.

50  : RichardsDensity(parameters),
51  _a(getParam<Real>("a")),
52  _b(getParam<Real>("b")),
53  _rt(getParam<Real>("temperature") * 8.314472), // multiply by gas constant
54  _molar_mass(getParam<Real>("molar_mass")),
55  _infinity_ratio(getParam<Real>("infinity_ratio")),
56  _rhs(_rt * _b / _a),
57  _b2oa(_b * _b / _a)
58 {
59  _vdw0 = densityVDW(0);
61 }
RichardsDensity(const InputParameters &parameters)
Real densityVDW(Real p) const
Density according to the van der Waals expression This is modified to yield density(p) as noted above...
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 _vdw0
density at P=0 according to the van der Waals expression
Real _rt
R*T (gas constant * temperature)
Real _molar_mass
molar mass of gas
Real ddensity(Real p) const
derivative of fluid density wrt porepressure

Member Function Documentation

◆ d2density()

Real RichardsDensityVDW::d2density ( Real  p) const
virtual

second derivative of fluid density wrt porepressure

Parameters
pporepressure

Implements RichardsDensity.

Definition at line 114 of file RichardsDensityVDW.C.

115 {
116  if (p >= 0)
117  {
118  Real y = _b2oa * p;
119  Real dy = _b2oa;
120  Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
121  Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
122  Real dsq = 0.5 / sq * (4.0 * 3.0 * Utility::pow<2>(-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy +
123  2.0 * (-2.0 - 18.0 * y + 9.0 * _rhs) * (-18.0 * dy));
124  Real d2sq = -dsq * dsq / sq;
125  d2sq += 0.5 / sq * (4.0 * 3.0 * 2.0 * (-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy * 3.0 * dy +
126  2.0 * (-18.0 * dy) * (-18.0 * dy));
127  Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
128  Real dcr =
129  1.0 / 3.0 * std::pow(-2.0 + 9.0 * _rhs - 18.0 * y + sq, -2.0 / 3.0) * (-18.0 * dy + dsq);
130  Real d2cr = 1.0 / 3.0 * (-2.0 / 3.0) * std::pow(-2 + 9 * _rhs - 18 * y + sq, -5. / 3.) *
131  Utility::pow<2>(-18 * dy + dsq);
132  d2cr += 1.0 / 3.0 * std::pow(-2 + 9 * _rhs - 18 * y + sq, -2. / 3.) * d2sq;
133  // Real dx = std::pow(2, 1.0/3.0)*( (3*dy)/3/cr + (-1 + 3*_rhs + 3*y)/3*(-dcr/cr/cr));
134  // dx -= dcr/3/std::pow(2, 1.0/3.0);
135  Real d2x = std::cbrt(2.0) *
136  (-(3.0 * dy) * dcr / (3.0 * cr * cr) + 3.0 * dy / 3.0 * (-dcr / (cr * cr)) +
137  (-1.0 + 3.0 * _rhs + 3.0 * y) / 3.0 *
138  (-d2cr / (cr * cr) + 2.0 * dcr * dcr / Utility::pow<3>(cr)));
139  d2x -= d2cr / (3.0 * std::cbrt(2.0));
140  return _molar_mass * d2x / _b;
141  }
142  else
143  return _infinity_ratio * _molar_mass * Utility::pow<2>(_slope0) * std::exp(_slope0 * p);
144 }
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

◆ ddensity()

Real RichardsDensityVDW::ddensity ( Real  p) const
virtual

derivative of fluid density wrt porepressure

Parameters
pporepressure

Implements RichardsDensity.

Definition at line 91 of file RichardsDensityVDW.C.

Referenced by RichardsDensityVDW().

92 {
93  if (p >= 0)
94  {
95  Real y = _b2oa * p;
96  Real dy = _b2oa;
97  Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
98  Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
99  Real dsq = 0.5 / sq * (4.0 * 3.0 * Utility::pow<2>(-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy +
100  2.0 * (-2.0 - 18.0 * y + 9.0 * _rhs) * (-18.0 * dy));
101  Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
102  Real dcr =
103  1.0 / 3.0 * std::pow(-2.0 + 9.0 * _rhs - 18.0 * y + sq, -2.0 / 3.0) * (-18.0 * dy + dsq);
104  Real dx = std::cbrt(2.0) *
105  ((3.0 * dy) / (3.0 * cr) + (-1.0 + 3.0 * _rhs + 3.0 * y) / 3.0 * (-dcr / (cr * cr)));
106  dx -= dcr / (3.0 * std::cbrt(2.0));
107  return _molar_mass * dx / _b;
108  }
109  else
110  return _infinity_ratio * _molar_mass * _slope0 * std::exp(_slope0 * p);
111 }
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

◆ density()

Real RichardsDensityVDW::density ( Real  p) const
virtual

fluid density as a function of porepressure

Parameters
pporepressure

Implements RichardsDensity.

Definition at line 82 of file RichardsDensityVDW.C.

83 {
84  if (p >= 0)
85  return densityVDW(p) - _vdw0;
86  else
87  return _infinity_ratio * _molar_mass * (std::exp(_slope0 * p) - 1.0);
88 }
Real densityVDW(Real p) const
Density according to the van der Waals expression This is modified to yield density(p) as noted above...
Real _infinity_ratio
density at P=-infinity is _infinity_ratio*_molar_mass
Real _slope0
(1/_molar_mass)*d(density)/dP at P=0
Real _vdw0
density at P=0 according to the van der Waals expression
Real _molar_mass
molar mass of gas

◆ densityVDW()

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 64 of file RichardsDensityVDW.C.

Referenced by density(), and RichardsDensityVDW().

65 {
66  // transform (P + a/V^2)(V-b) = RT to
67  // (y + x^2)(1/x - 1) = rhs
68  // with: y = b^2*P/a, x = b/V, rhs = RT*b/a
69  // Then solve for x, using mathematica
70  // Then density = molar_mass/V = molar_mass*x/b
71  Real y = _b2oa * p;
72  Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
73  Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
74  Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
75  Real x = 1.0 / 3.0;
76  x += std::cbrt(2.0) * (-1.0 + 3.0 * _rhs + 3.0 * y) / (3.0 * cr);
77  x -= cr / (3.0 * std::cbrt(2.0));
78  return _molar_mass * x / _b;
79 }
Real _b
van der Waals b
Real _molar_mass
molar mass of gas

◆ execute()

void RichardsDensity::execute ( )
inherited

Definition at line 34 of file RichardsDensity.C.

35 {
36 }

◆ finalize()

void RichardsDensity::finalize ( )
inherited

Definition at line 39 of file RichardsDensity.C.

40 {
41 }

◆ initialize()

void RichardsDensity::initialize ( )
inherited

Definition at line 29 of file RichardsDensity.C.

30 {
31 }

Member Data Documentation

◆ _a

Real RichardsDensityVDW::_a
protected

van der Waals a

Definition at line 62 of file RichardsDensityVDW.h.

◆ _b

Real RichardsDensityVDW::_b
protected

van der Waals b

Definition at line 65 of file RichardsDensityVDW.h.

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

◆ _b2oa

Real RichardsDensityVDW::_b2oa
protected

b^2/a

Definition at line 80 of file RichardsDensityVDW.h.

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

◆ _infinity_ratio

Real RichardsDensityVDW::_infinity_ratio
protected

density at P=-infinity is _infinity_ratio*_molar_mass

Definition at line 74 of file RichardsDensityVDW.h.

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

◆ _molar_mass

Real RichardsDensityVDW::_molar_mass
protected

molar mass of gas

Definition at line 71 of file RichardsDensityVDW.h.

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

◆ _rhs

Real RichardsDensityVDW::_rhs
protected

R*T*b/a.

Definition at line 77 of file RichardsDensityVDW.h.

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

◆ _rt

Real RichardsDensityVDW::_rt
protected

R*T (gas constant * temperature)

Definition at line 68 of file RichardsDensityVDW.h.

◆ _slope0

Real RichardsDensityVDW::_slope0
protected

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

Definition at line 86 of file RichardsDensityVDW.h.

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

◆ _vdw0

Real RichardsDensityVDW::_vdw0
protected

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

Definition at line 83 of file RichardsDensityVDW.h.

Referenced by density(), and RichardsDensityVDW().


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