www.mooseframework.org
PorousFlowCapillaryPressure.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
11 
14 {
16  params.addRangeCheckedParam<Real>(
17  "sat_lr",
18  0.0,
19  "sat_lr >= 0 & sat_lr < 1",
20  "Liquid residual saturation. Must be between 0 and 1. Default is 0");
21  params.addRangeCheckedParam<Real>("pc_max",
22  1.0e9,
23  "pc_max > 0",
24  "Maximum capillary pressure (Pa). Must be > 0. Default is 1e9");
25  params.addParam<bool>("log_extension",
26  true,
27  "Use a logarithmic extension for low saturation to avoid capillary "
28  "pressure going to infinity. Default is true. Set to false if your "
29  "capillary pressure depends on spatially-dependent variables other than "
30  "saturation, as the log-extension C++ code for this case has yet to be "
31  "implemented");
32  params.addClassDescription("Capillary pressure base class");
33  return params;
34 }
35 
37  : DiscreteElementUserObject(parameters),
38  _sat_lr(getParam<Real>("sat_lr")),
39  _dseff_ds(1.0 / (1.0 - _sat_lr)),
40  _log_ext(getParam<bool>("log_extension")),
41  _pc_max(getParam<Real>("pc_max")),
42  _sat_ext(0.0),
43  _pc_ext(0.0),
44  _slope_ext(0.0),
45  _log10(std::log(10.0))
46 {
47 }
48 
49 void
51 {
52  // If _log_ext = true, calculate the saturation where the the logarithmic
53  // extension meets the raw capillary pressure curve.
54  if (_log_ext)
55  {
58  _slope_ext = (std::log10(_pc_ext) - std::log10(_pc_max)) / _sat_ext;
59  }
60 }
61 
62 Real
63 PorousFlowCapillaryPressure::capillaryPressure(Real saturation, unsigned qp) const
64 {
65  if (_log_ext && saturation < _sat_ext)
67  else
69 }
70 
71 ADReal
72 PorousFlowCapillaryPressure::capillaryPressure(const ADReal & saturation, unsigned qp) const
73 {
74  const Real Pc = capillaryPressure(saturation.value(), qp);
75  const Real dPc_ds = dCapillaryPressure(saturation.value(), qp);
76 
77  ADReal result = Pc;
78  result.derivatives() = saturation.derivatives() * dPc_ds;
79 
80  return result;
81 }
82 
83 Real
84 PorousFlowCapillaryPressure::dCapillaryPressure(Real saturation, unsigned qp) const
85 {
86  if (_log_ext && saturation < _sat_ext)
88  else
90 }
91 
92 ADReal
93 PorousFlowCapillaryPressure::dCapillaryPressure(const ADReal & saturation, unsigned qp) const
94 {
95  const Real dPc = dCapillaryPressure(saturation.value(), qp);
96  const Real d2Pc_ds2 = d2CapillaryPressure(saturation.value(), qp);
97 
98  ADReal result = dPc;
99  result.derivatives() = saturation.derivatives() * d2Pc_ds2;
100 
101  return result;
102 }
103 
104 Real
105 PorousFlowCapillaryPressure::d2CapillaryPressure(Real saturation, unsigned qp) const
106 {
107  if (_log_ext && saturation < _sat_ext)
109  else
111 }
112 
113 Real
115 {
116  return (saturation - _sat_lr) / (1.0 - _sat_lr);
117 }
118 
119 Real
120 PorousFlowCapillaryPressure::saturation(Real pc, unsigned qp) const
121 {
122  return effectiveSaturation(pc, qp) * (1.0 - _sat_lr) + _sat_lr;
123 }
124 
125 ADReal
126 PorousFlowCapillaryPressure::saturation(const ADReal & pc, unsigned qp) const
127 {
128  const Real s = effectiveSaturation(pc.value(), qp) * (1.0 - _sat_lr) + _sat_lr;
129  const Real ds_dpc = dEffectiveSaturation(pc.value(), qp) * (1.0 - _sat_lr);
130 
131  ADReal result = s;
132  result.derivatives() = pc.derivatives() * ds_dpc;
133 
134  return result;
135 }
136 
137 Real
138 PorousFlowCapillaryPressure::dSaturation(Real pc, unsigned qp) const
139 {
140  return dEffectiveSaturation(pc, qp) * (1.0 - _sat_lr);
141 }
142 
143 ADReal
144 PorousFlowCapillaryPressure::dSaturation(const ADReal & pc, unsigned qp) const
145 {
146  const Real ds = dEffectiveSaturation(pc.value(), qp) * (1.0 - _sat_lr);
147  const Real d2s_dpc2 = d2EffectiveSaturation(pc.value(), qp) * (1.0 - _sat_lr);
148 
149  ADReal result = ds;
150  result.derivatives() = pc.derivatives() * d2s_dpc2;
151 
152  return result;
153 }
154 
155 Real
156 PorousFlowCapillaryPressure::d2Saturation(Real pc, unsigned qp) const
157 {
158  return d2EffectiveSaturation(pc, qp) * (1.0 - _sat_lr);
159 }
160 
161 Real
163 {
164  return _pc_ext * std::pow(10.0, _slope_ext * (saturation - _sat_ext));
165 }
166 
167 Real
169 {
170  return _pc_ext * _slope_ext * _log10 * std::pow(10.0, _slope_ext * (saturation - _sat_ext));
171 }
172 
173 Real
175 {
176  return _pc_ext * _slope_ext * _slope_ext * _log10 * _log10 *
177  std::pow(10.0, _slope_ext * (saturation - _sat_ext));
178 }
179 
180 Real
182 {
183  // Initial guess for saturation where extension matches curve
184  Real sat = _sat_lr + 0.01;
185 
186  // Calculate the saturation where the extension matches the derivative of the
187  // raw capillary pressure curve
188  const unsigned int max_its = 20;
189  const Real nr_tol = 1.0e-8;
190  unsigned int iter = 0;
191 
192  while (std::abs(interceptFunction(sat)) > nr_tol)
193  {
194  sat = sat - interceptFunction(sat) / interceptFunctionDeriv(sat);
195 
196  iter++;
197  if (iter > max_its)
198  break;
199  }
200 
201  return sat;
202 }
203 
204 Real
206 {
209 
210  return std::log10(pc) - saturation * dpc / (_log10 * pc) - std::log10(_pc_max);
211 }
212 
213 Real
215 {
219 
220  return saturation * (dpc * dpc / pc - d2pc) / (_log10 * pc);
221 }
virtual Real dCapillaryPressureCurve(Real saturation, unsigned qp=0) const =0
Derivative of raw capillary pressure wrt true saturation.
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
virtual Real effectiveSaturation(Real pc, unsigned qp=0) const =0
Effective saturation as a function of capillary pressure.
virtual Real d2EffectiveSaturation(Real pc, unsigned qp=0) const =0
Second derivative of effective saturation wrt capillary pressure.
static InputParameters validParams()
Real effectiveSaturationFromSaturation(Real saturation) const
Effective saturation of liquid phase given liquid saturation and residual liquid saturation.
Real saturation(Real pc, unsigned qp=0) const
Saturation as a function of capillary pressure.
virtual Real capillaryPressureCurve(Real saturation, unsigned qp=0) const =0
Raw capillary pressure curve (does not include logarithmic extension)
Real _slope_ext
Gradient of the logarithmic extension This is computed only for qp=0.
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
const Real _pc_max
Maximum capillary pressure (Pa). Note: must be <= 0.
Real interceptFunctionDeriv(Real s) const
Calculates the saturation where the logarithmic extension to capillary pressure at low saturation...
Real extensionSaturation() const
Calculates the saturation where the logarithmic extension to capillary pressure meets the raw curve u...
virtual Real capillaryPressure(Real saturation, unsigned qp=0) const
Capillary pressure is calculated as a function of true saturation.
Real dSaturation(Real pc, unsigned qp=0) const
Derivative of saturation wrt capillary pressure.
DualReal ADReal
auto log(const T &)
Real d2CapillaryPressureLogExt(Real s) const
The second derivative of capillary pressure in the logarithmic extension This implementation assumes ...
PorousFlowCapillaryPressure(const InputParameters &parameters)
Real _pc_ext
Capillary pressure where the extension meets the raw curve.
virtual Real d2CapillaryPressure(Real saturation, unsigned qp=0) const
Second derivative of capillary pressure wrt true saturation.
virtual Real d2CapillaryPressureCurve(Real saturation, unsigned qp=0) const =0
Second derivative of raw capillary pressure wrt true saturation.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real interceptFunction(Real s) const
Calculates the saturation where the logarithmic extension to capillary pressure at low saturation...
Real capillaryPressureLogExt(Real s) const
The capillary pressure in the logarithmic extension This implementation assumes capillary-pressure is...
virtual Real dCapillaryPressure(Real saturation, unsigned qp=0) const
Derivative of capillary pressure wrt true saturation.
Real _sat_ext
Saturation where the logarithmic extension meets the raw curve This is computed only for qp=0...
const Real _sat_lr
Liquid residual saturation.
void addClassDescription(const std::string &doc_string)
bool _log_ext
Flag to use a logarithmic extension for low saturation.
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
Real d2Saturation(Real pc, unsigned qp=0) const
Second derivative of saturation wrt capillary pressure.
static InputParameters validParams()
MooseUnits pow(const MooseUnits &, int)
Real dCapillaryPressureLogExt(Real s) const
The derivative of capillary pressure in the logarithmic extension This implementation assumes capilla...
virtual Real dEffectiveSaturation(Real pc, unsigned qp=0) const =0
Derivative of effective saturation wrt capillary pressure.