www.mooseframework.org
PorousFlowCapillaryPressure.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
9 
10 template <>
11 InputParameters
13 {
14  InputParameters params = validParams<GeneralUserObject>();
15  params.addRangeCheckedParam<Real>(
16  "sat_lr",
17  0.0,
18  "sat_lr >= 0 & sat_lr <= 1",
19  "Liquid residual saturation. Must be between 0 and 1. Default is 0");
20  params.addRangeCheckedParam<Real>("pc_max",
21  1.0e9,
22  "pc_max > 0",
23  "Maximum capillary pressure (Pa). Must be > 0. Default is 1e9");
24  params.addParam<bool>("log_extension",
25  true,
26  "Use a logarithmic extension for low saturation to avoid capillary "
27  "pressure going to infinity. Default is true");
28  params.addClassDescription("Capillary pressure base class");
29  return params;
30 }
31 
33  : GeneralUserObject(parameters),
34  _sat_lr(getParam<Real>("sat_lr")),
35  _dseff_ds(1.0 / (1.0 - _sat_lr)),
36  _log_ext(getParam<bool>("log_extension")),
37  _pc_max(getParam<Real>("pc_max")),
38  _sat_ext(0.0),
39  _pc_ext(0.0),
40  _slope_ext(0.0),
41  _log10(std::log(10.0))
42 {
43 }
44 
45 void
47 {
48  // If _log_ext = true, calculate the saturation where the the logarithmic
49  // extension meets the raw capillary pressure curve.
50  if (_log_ext)
51  {
54  _slope_ext = (std::log10(_pc_ext) - std::log10(_pc_max)) / _sat_ext;
55  }
56 }
57 
58 Real
60 {
61  if (_log_ext && saturation < _sat_ext)
62  return capillaryPressureLogExt(saturation);
63  else
64  return capillaryPressureCurve(saturation);
65 }
66 
67 Real
69 {
70  if (_log_ext && saturation < _sat_ext)
71  return dCapillaryPressureLogExt(saturation);
72  else
73  return dCapillaryPressureCurve(saturation);
74 }
75 
76 Real
78 {
79  if (_log_ext && saturation < _sat_ext)
80  return d2CapillaryPressureLogExt(saturation);
81  else
82  return d2CapillaryPressureCurve(saturation);
83 }
84 
85 Real
87 {
88  return (saturation - _sat_lr) / (1.0 - _sat_lr);
89 }
90 
91 Real
93 {
94  return _pc_ext * std::pow(10.0, _slope_ext * (saturation - _sat_ext));
95 }
96 
97 Real
99 {
100  return _pc_ext * _slope_ext * _log10 * std::pow(10.0, _slope_ext * (saturation - _sat_ext));
101 }
102 
103 Real
105 {
106  return _pc_ext * _slope_ext * _slope_ext * _log10 * _log10 *
107  std::pow(10.0, _slope_ext * (saturation - _sat_ext));
108 }
109 
110 Real
112 {
113  // Initial guess for saturation where extension matches curve
114  Real sat = _sat_lr + 0.01;
115 
116  // Calculate the saturation where the extension matches the derivative of the
117  // raw capillary pressure curve
118  const unsigned int max_its = 20;
119  const Real nr_tol = 1.0e-8;
120  unsigned int iter = 0;
121 
122  while (std::abs(interceptFunction(sat)) > nr_tol)
123  {
124  sat = sat - interceptFunction(sat) / interceptFunctionDeriv(sat);
125 
126  iter++;
127  if (iter > max_its)
128  break;
129  }
130 
131  return sat;
132 }
133 
134 Real
136 {
137  Real pc = capillaryPressureCurve(saturation);
138  Real dpc = dCapillaryPressureCurve(saturation);
139 
140  return std::log10(pc) - saturation * dpc / (_log10 * pc) - std::log10(_pc_max);
141 }
142 
143 Real
145 {
146  Real pc = capillaryPressureCurve(saturation);
147  Real dpc = dCapillaryPressureCurve(saturation);
148  Real d2pc = d2CapillaryPressureCurve(saturation);
149 
150  return saturation * (dpc * dpc / pc - d2pc) / (_log10 * pc);
151 }
Real d2CapillaryPressureLogExt(Real s) const
The second derivative of capillary pressure in the logarithmic extension.
virtual Real dCapillaryPressureCurve(Real saturation) const =0
Derivative of raw capillary pressure wrt true saturation.
virtual Real dCapillaryPressure(Real saturation) const
Derivative of capillary pressure wrt true saturation.
Real _slope_ext
Gradient of the logarithmic extension.
Real extensionSaturation() const
Calculates the saturation where the logarithmic extension to capillary pressure meets the raw curve u...
InputParameters validParams< PorousFlowCapillaryPressure >()
virtual Real d2CapillaryPressure(Real saturation) const
Second derivative of capillary pressure wrt true saturation.
const Real _pc_max
Maximum capillary pressure (Pa). Note: must be <= 0.
virtual Real capillaryPressure(Real saturation) const
Capillary pressure is calculated as a function of true saturation.
Real capillaryPressureLogExt(Real s) const
The capillary pressure in the logarithmic extension.
Real effectiveSaturationFromSaturation(Real saturation) const
Effective saturation of liquid phase given liquid saturation and residual liquid saturation.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
PorousFlowCapillaryPressure(const InputParameters &parameters)
virtual Real d2CapillaryPressureCurve(Real saturation) const =0
Second derivative of raw capillary pressure wrt true saturation.
Real _pc_ext
Capillary pressure where the extension meets the raw curve.
Real interceptFunction(Real s) const
Calculates the saturation where the logarithmic extension to capillary pressure at low saturation...
Real _sat_ext
Saturation where the logarithmic extension meets the raw curve.
const Real _sat_lr
Liquid residual saturation.
Real interceptFunctionDeriv(Real s) const
Calculates the saturation where the logarithmic extension to capillary pressure at low saturation...
bool _log_ext
Flag to use a logarithmic extension for low saturation.
virtual Real capillaryPressureCurve(Real saturation) const =0
Raw capillary pressure curve (does not include logarithmic extension)
void FORTRAN_CALL() saturation(double &P, double &T, int &N, int &nerr)
Real dCapillaryPressureLogExt(Real s) const
The derivative of capillary pressure in the logarithmic extension.