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             : //  Fluid density assuming constant bulk modulus, but reducing via a cubic to zero at p=zero_point
      11             : //
      12             : #include "RichardsDensityConstBulkCut.h"
      13             : 
      14             : registerMooseObject("RichardsApp", RichardsDensityConstBulkCut);
      15             : 
      16             : InputParameters
      17           7 : RichardsDensityConstBulkCut::validParams()
      18             : {
      19           7 :   InputParameters params = RichardsDensity::validParams();
      20          14 :   params.addRequiredParam<Real>("dens0", "Reference density of fluid.  Eg 1000");
      21          14 :   params.addRequiredParam<Real>("bulk_mod", "Bulk modulus of fluid.  Eg 2E9");
      22          14 :   params.addParam<Real>("cut_limit",
      23          14 :                         1E5,
      24             :                         "For porepressure > cut_limit, density = dens0*Exp(pressure/bulk).  For "
      25             :                         "porepressure < cut_limie, density = cubic*dens0*Exp(pressure/bulk), where "
      26             :                         "cubic=1 for pressure=cut_limit, and cubic=0 for pressure<=zero_point");
      27          14 :   params.addParam<Real>("zero_point",
      28          14 :                         0,
      29             :                         "For porepressure > cut_limit, density = dens0*Exp(pressure/bulk).  For "
      30             :                         "porepressure < cut_limie, density = cubic*dens0*Exp(pressure/bulk), where "
      31             :                         "cubic=1 for pressure=cut_limit, and cubic=0 for pressure<=zero_point");
      32           7 :   params.addClassDescription(
      33             :       "Fluid density assuming constant bulk modulus.  dens0 * Exp(pressure/bulk)");
      34           7 :   return params;
      35           0 : }
      36             : 
      37           4 : RichardsDensityConstBulkCut::RichardsDensityConstBulkCut(const InputParameters & parameters)
      38             :   : RichardsDensity(parameters),
      39           4 :     _dens0(getParam<Real>("dens0")),
      40           8 :     _bulk(getParam<Real>("bulk_mod")),
      41           8 :     _cut_limit(getParam<Real>("cut_limit")),
      42           8 :     _zero_point(getParam<Real>("zero_point")),
      43           8 :     _c3(std::pow(_cut_limit - _zero_point, 3))
      44             : {
      45           4 :   if (_zero_point >= _cut_limit)
      46           1 :     mooseError("RichardsDensityConstantbulkCut: zero_point must be less than cut_limit");
      47           3 : }
      48             : 
      49             : Real
      50        9230 : RichardsDensityConstBulkCut::density(Real p) const
      51             : {
      52        9230 :   if (p <= _zero_point)
      53             :     return 0;
      54             : 
      55        8884 :   Real unmodified = _dens0 * std::exp(p / _bulk);
      56        8884 :   if (p >= _cut_limit)
      57             :     return unmodified;
      58             : 
      59        8538 :   Real modifier =
      60        8538 :       (3 * _cut_limit - 2 * p - _zero_point) * (p - _zero_point) * (p - _zero_point) / _c3;
      61        8538 :   return modifier * unmodified;
      62             : }
      63             : 
      64             : Real
      65        1934 : RichardsDensityConstBulkCut::ddensity(Real p) const
      66             : {
      67        1934 :   if (p <= _zero_point)
      68             :     return 0;
      69             : 
      70        1684 :   Real unmodified = _dens0 * std::exp(p / _bulk);
      71        1684 :   if (p >= _cut_limit)
      72         250 :     return unmodified / _bulk;
      73             : 
      74        1434 :   Real modifier =
      75        1434 :       (3 * _cut_limit - 2 * p - _zero_point) * (p - _zero_point) * (p - _zero_point) / _c3;
      76        1434 :   Real dmodifier = 6 * (_cut_limit - p) * (p - _zero_point) / _c3;
      77        1434 :   return (modifier / _bulk + dmodifier) * unmodified;
      78             : }
      79             : 
      80             : Real
      81         606 : RichardsDensityConstBulkCut::d2density(Real p) const
      82             : {
      83         606 :   if (p <= _zero_point)
      84             :     return 0;
      85             : 
      86         360 :   Real unmodified = _dens0 * std::exp(p / _bulk);
      87         360 :   if (p >= _cut_limit)
      88         246 :     return unmodified / _bulk / _bulk;
      89             : 
      90         114 :   Real modifier =
      91         114 :       (3 * _cut_limit - 2 * p - _zero_point) * (p - _zero_point) * (p - _zero_point) / _c3;
      92         114 :   Real dmodifier = 6 * (_cut_limit - p) * (p - _zero_point) / _c3;
      93         114 :   Real d2modifier = (6 * _cut_limit + 6 * _zero_point - 12 * p) / _c3;
      94         114 :   return (modifier / _bulk / _bulk + 2 * dmodifier / _bulk + d2modifier) * unmodified;
      95             : }