Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "RichardsDensityVDW.h"
      11             : #include "libmesh/utility.h"
      12             : 
      13             : registerMooseObject("RichardsApp", RichardsDensityVDW);
      14             : 
      15             : InputParameters
      16           8 : RichardsDensityVDW::validParams()
      17             : {
      18           8 :   InputParameters params = RichardsDensity::validParams();
      19          16 :   params.addRequiredRangeCheckedParam<Real>(
      20             :       "a",
      21             :       "a > 0",
      22             :       "Parameter 'a' in the van der Waals expression (P + n^2 a/V^2)(V - nb) = nRT.  "
      23             :       "Example for methane 0.2303 Pa m^6 mol^-2");
      24          16 :   params.addRequiredRangeCheckedParam<Real>(
      25             :       "b",
      26             :       "b > 0",
      27             :       "Parameter 'b' in the van der Waals expression (P + n^2 a/V^2)(V - nb) = nRT.  "
      28             :       "Example for methane 4.31E-5 m^3/mol");
      29          16 :   params.addRequiredRangeCheckedParam<Real>(
      30             :       "temperature", "temperature > 0", "Temperature in Kelvin");
      31          16 :   params.addRequiredRangeCheckedParam<Real>(
      32             :       "molar_mass",
      33             :       "molar_mass > 0",
      34             :       "Molar mass of the gas.  Example for methane 16.04246E-3 kg/mol");
      35          16 :   params.addRangeCheckedParam<Real>("infinity_ratio",
      36             :                                     10,
      37             :                                     "infinity_ratio > 0",
      38             :                                     "For P<0 the density is not physically defined, "
      39             :                                     "but numerically it is advantageous to define "
      40             :                                     "it:  density(P=-infinity) = "
      41             :                                     "-infinity_ratio*molar_mass, and density tends "
      42             :                                     "exponentially towards this value as P -> "
      43             :                                     "-infinity.  (Units are mol/m^3).");
      44           8 :   params.addClassDescription("Density of van der Waals gas.");
      45           8 :   return params;
      46           0 : }
      47             : 
      48           4 : RichardsDensityVDW::RichardsDensityVDW(const InputParameters & parameters)
      49             :   : RichardsDensity(parameters),
      50           4 :     _a(getParam<Real>("a")),
      51           8 :     _b(getParam<Real>("b")),
      52           8 :     _rt(getParam<Real>("temperature") * 8.314472), // multiply by gas constant
      53           8 :     _molar_mass(getParam<Real>("molar_mass")),
      54           8 :     _infinity_ratio(getParam<Real>("infinity_ratio")),
      55           4 :     _rhs(_rt * _b / _a),
      56           4 :     _b2oa(_b * _b / _a)
      57             : {
      58           4 :   _vdw0 = densityVDW(0);
      59           4 :   _slope0 = ddensity(0) / (_molar_mass * _infinity_ratio);
      60           4 : }
      61             : 
      62             : Real
      63         310 : RichardsDensityVDW::densityVDW(Real p) const
      64             : {
      65             :   // transform (P + a/V^2)(V-b) = RT to
      66             :   // (y + x^2)(1/x - 1) = rhs
      67             :   // with: y = b^2*P/a, x = b/V, rhs = RT*b/a
      68             :   // Then solve for x, using mathematica
      69             :   // Then density = molar_mass/V = molar_mass*x/b
      70         310 :   Real y = _b2oa * p;
      71         310 :   Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
      72         310 :                       Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
      73         310 :   Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
      74             :   Real x = 1.0 / 3.0;
      75         310 :   x += std::cbrt(2.0) * (-1.0 + 3.0 * _rhs + 3.0 * y) / (3.0 * cr);
      76         310 :   x -= cr / (3.0 * std::cbrt(2.0));
      77         310 :   return _molar_mass * x / _b;
      78             : }
      79             : 
      80             : Real
      81         606 : RichardsDensityVDW::density(Real p) const
      82             : {
      83         606 :   if (p >= 0)
      84         306 :     return densityVDW(p) - _vdw0;
      85             :   else
      86         300 :     return _infinity_ratio * _molar_mass * (std::exp(_slope0 * p) - 1.0);
      87             : }
      88             : 
      89             : Real
      90         610 : RichardsDensityVDW::ddensity(Real p) const
      91             : {
      92         610 :   if (p >= 0)
      93             :   {
      94         310 :     Real y = _b2oa * p;
      95             :     Real dy = _b2oa;
      96         310 :     Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
      97         310 :                         Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
      98         310 :     Real dsq = 0.5 / sq *
      99         310 :                (4.0 * 3.0 * Utility::pow<2>(-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy +
     100         310 :                 2.0 * (-2.0 - 18.0 * y + 9.0 * _rhs) * (-18.0 * dy));
     101         310 :     Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
     102             :     Real dcr =
     103         310 :         1.0 / 3.0 * std::pow(-2.0 + 9.0 * _rhs - 18.0 * y + sq, -2.0 / 3.0) * (-18.0 * dy + dsq);
     104         310 :     Real dx = std::cbrt(2.0) *
     105         310 :               ((3.0 * dy) / (3.0 * cr) + (-1.0 + 3.0 * _rhs + 3.0 * y) / 3.0 * (-dcr / (cr * cr)));
     106         310 :     dx -= dcr / (3.0 * std::cbrt(2.0));
     107         310 :     return _molar_mass * dx / _b;
     108             :   }
     109             :   else
     110         300 :     return _infinity_ratio * _molar_mass * _slope0 * std::exp(_slope0 * p);
     111             : }
     112             : 
     113             : Real
     114         606 : RichardsDensityVDW::d2density(Real p) const
     115             : {
     116         606 :   if (p >= 0)
     117             :   {
     118         306 :     Real y = _b2oa * p;
     119             :     Real dy = _b2oa;
     120         306 :     Real sq = std::sqrt(4.0 * Utility::pow<3>(-1.0 + 3.0 * y + 3.0 * _rhs) +
     121         306 :                         Utility::pow<2>(-2.0 - 18.0 * y + 9.0 * _rhs));
     122         306 :     Real dsq = 0.5 / sq *
     123         306 :                (4.0 * 3.0 * Utility::pow<2>(-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy +
     124         306 :                 2.0 * (-2.0 - 18.0 * y + 9.0 * _rhs) * (-18.0 * dy));
     125         306 :     Real d2sq = -dsq * dsq / sq;
     126         306 :     d2sq += 0.5 / sq *
     127         306 :             (4.0 * 3.0 * 2.0 * (-1.0 + 3.0 * y + 3.0 * _rhs) * 3.0 * dy * 3.0 * dy +
     128         306 :              2.0 * (-18.0 * dy) * (-18.0 * dy));
     129         306 :     Real cr = std::cbrt(-2.0 + 9.0 * _rhs - 18.0 * y + sq);
     130             :     Real dcr =
     131         306 :         1.0 / 3.0 * std::pow(-2.0 + 9.0 * _rhs - 18.0 * y + sq, -2.0 / 3.0) * (-18.0 * dy + dsq);
     132         306 :     Real d2cr = 1.0 / 3.0 * (-2.0 / 3.0) * std::pow(-2 + 9 * _rhs - 18 * y + sq, -5. / 3.) *
     133         306 :                 Utility::pow<2>(-18 * dy + dsq);
     134         306 :     d2cr += 1.0 / 3.0 * std::pow(-2 + 9 * _rhs - 18 * y + sq, -2. / 3.) * d2sq;
     135             :     // Real dx = std::pow(2, 1.0/3.0)*( (3*dy)/3/cr + (-1 + 3*_rhs + 3*y)/3*(-dcr/cr/cr));
     136             :     // dx -= dcr/3/std::pow(2, 1.0/3.0);
     137             :     Real d2x = std::cbrt(2.0) *
     138         306 :                (-(3.0 * dy) * dcr / (3.0 * cr * cr) + 3.0 * dy / 3.0 * (-dcr / (cr * cr)) +
     139         306 :                 (-1.0 + 3.0 * _rhs + 3.0 * y) / 3.0 *
     140         306 :                     (-d2cr / (cr * cr) + 2.0 * dcr * dcr / Utility::pow<3>(cr)));
     141         306 :     d2x -= d2cr / (3.0 * std::cbrt(2.0));
     142         306 :     return _molar_mass * d2x / _b;
     143             :   }
     144             :   else
     145         300 :     return _infinity_ratio * _molar_mass * Utility::pow<2>(_slope0) * std::exp(_slope0 * p);
     146             : }