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
	RT (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
	RTb/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 36 of file RichardsDensityVDW.h.

Constructor & Destructor Documentation

◆ RichardsDensityVDW()

RichardsDensityVDW::RichardsDensityVDW ( const InputParameters & parameters )

Definition at line 49 of file RichardsDensityVDW.C.

   : RichardsDensity(parameters),
     _a(getParam<Real>("a")),
     _b(getParam<Real>("b")),
     _rt(getParam<Real>("temperature") * 8.314472), // multiply by gas constant
     _molar_mass(getParam<Real>("molar_mass")),
     _infinity_ratio(getParam<Real>("infinity_ratio")),
     _rhs(_rt * _b / _a),
     _b2oa(_b * _b / _a)
 {
   _vdw0 = densityVDW(0);
   _slope0 = ddensity(0) / (_molar_mass * _infinity_ratio);
 }

Member Function Documentation

◆ d2density()

Real RichardsDensityVDW::d2density ( Real p ) const

virtual

second derivative of fluid density wrt porepressure

Parameters

p	porepressure

Implements RichardsDensity.

Definition at line 114 of file RichardsDensityVDW.C.

 {
   if (p >= 0)
   {
     Real y = _b2oa * p;
     Real dy = _b2oa;
     Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
                         Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
     Real dsq = 0.5 / sq * (4.0 * 3.0 * Utility::pow<2>(-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy +
                            2.0 * (-2.0 - 18.0 * y + 9.0 * _rhs) * (-18.0 * dy));
     Real d2sq = -dsq * dsq / sq;
     d2sq += 0.5 / sq * (4.0 * 3.0 * 2.0 * (-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy * 3.0 * dy +
                         2.0 * (-18.0 * dy) * (-18.0 * dy));
     Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
     Real dcr =
         1.0 / 3.0 * std::pow(-2.0 + 9.0 * _rhs - 18.0 * y + sq, -2.0 / 3.0) * (-18.0 * dy + dsq);
     Real d2cr = 1.0 / 3.0 * (-2.0 / 3.0) * std::pow(-2 + 9 * _rhs - 18 * y + sq, -5. / 3.) *
                 Utility::pow<2>(-18 * dy + dsq);
     d2cr += 1.0 / 3.0 * std::pow(-2 + 9 * _rhs - 18 * y + sq, -2. / 3.) * d2sq;
     // Real dx = std::pow(2, 1.0/3.0)*( (3*dy)/3/cr + (-1 + 3*_rhs + 3*y)/3*(-dcr/cr/cr));
     // dx -= dcr/3/std::pow(2, 1.0/3.0);
     Real d2x = std::cbrt(2.0) *
                (-(3.0 * dy) * dcr / (3.0 * cr * cr) + 3.0 * dy / 3.0 * (-dcr / (cr * cr)) +
                 (-1.0 + 3.0 * _rhs + 3.0 * y) / 3.0 *
                     (-d2cr / (cr * cr) + 2.0 * dcr * dcr / Utility::pow<3>(cr)));
     d2x -= d2cr / (3.0 * std::cbrt(2.0));
     return _molar_mass * d2x / _b;
   }
   else
     return _infinity_ratio * _molar_mass * Utility::pow<2>(_slope0) * std::exp(_slope0 * p);
 }

◆ ddensity()

Real RichardsDensityVDW::ddensity ( Real p ) const

virtual

derivative of fluid density wrt porepressure

Parameters

p	porepressure

Implements RichardsDensity.

Definition at line 91 of file RichardsDensityVDW.C.

 {
   if (p >= 0)
   {
     Real y = _b2oa * p;
     Real dy = _b2oa;
     Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
                         Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
     Real dsq = 0.5 / sq * (4.0 * 3.0 * Utility::pow<2>(-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy +
                            2.0 * (-2.0 - 18.0 * y + 9.0 * _rhs) * (-18.0 * dy));
     Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
     Real dcr =
         1.0 / 3.0 * std::pow(-2.0 + 9.0 * _rhs - 18.0 * y + sq, -2.0 / 3.0) * (-18.0 * dy + dsq);
     Real dx = std::cbrt(2.0) *
               ((3.0 * dy) / (3.0 * cr) + (-1.0 + 3.0 * _rhs + 3.0 * y) / 3.0 * (-dcr / (cr * cr)));
     dx -= dcr / (3.0 * std::cbrt(2.0));
     return _molar_mass * dx / _b;
   }
   else
     return _infinity_ratio * _molar_mass * _slope0 * std::exp(_slope0 * p);
 }

Referenced by RichardsDensityVDW().

◆ density()

Real RichardsDensityVDW::density ( Real p ) const

virtual

fluid density as a function of porepressure

Parameters

p	porepressure

Implements RichardsDensity.

Definition at line 82 of file RichardsDensityVDW.C.

 {
   if (p >= 0)
     return densityVDW(p) - _vdw0;
   else
     return _infinity_ratio * _molar_mass * (std::exp(_slope0 * p) - 1.0);
 }

◆ 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

p	porepressure

Definition at line 64 of file RichardsDensityVDW.C.

 {
   // transform (P + a/V^2)(V-b) = RT to
   // (y + x^2)(1/x - 1) = rhs
   // with: y = b^2*P/a, x = b/V, rhs = RT*b/a
   // Then solve for x, using mathematica
   // Then density = molar_mass/V = molar_mass*x/b
   Real y = _b2oa * p;
   Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
                       Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
   Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
   Real x = 1.0 / 3.0;
   x += std::cbrt(2.0) * (-1.0 + 3.0 * _rhs + 3.0 * y) / (3.0 * cr);
   x -= cr / (3.0 * std::cbrt(2.0));
   return _molar_mass * x / _b;
 }

Referenced by density(), and RichardsDensityVDW().

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

◆ _b

Real RichardsDensityVDW::_b

protected

van der Waals b

Definition at line 64 of file RichardsDensityVDW.h.

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

◆ _b2oa

Real RichardsDensityVDW::_b2oa

protected

b^2/a

Definition at line 79 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 73 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 70 of file RichardsDensityVDW.h.

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

◆ _rhs

Real RichardsDensityVDW::_rhs

protected

R*T*b/a.

Definition at line 76 of file RichardsDensityVDW.h.

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

◆ _rt

Real RichardsDensityVDW::_rt

protected

R*T (gas constant * temperature)

Definition at line 67 of file RichardsDensityVDW.h.

◆ _slope0

Real RichardsDensityVDW::_slope0

protected

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

Definition at line 85 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 82 of file RichardsDensityVDW.h.

Referenced by density(), and RichardsDensityVDW().

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

richards/include/userobjects/RichardsDensityVDW.h
richards/src/userobjects/RichardsDensityVDW.C

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ RichardsDensityVDW()

Member Function Documentation

◆ d2density()

◆ ddensity()

◆ density()

◆ densityVDW()

◆ execute()

◆ finalize()

◆ initialize()

Member Data Documentation

◆ _a

◆ _b

◆ _b2oa

◆ _infinity_ratio

◆ _molar_mass

◆ _rhs

◆ _rt

◆ _slope0

◆ _vdw0