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 "PorousFlowCapillaryPressure.h"
11 :
12 : InputParameters
13 11411 : PorousFlowCapillaryPressure::validParams()
14 : {
15 11411 : InputParameters params = DiscreteElementUserObject::validParams();
16 34233 : params.addRangeCheckedParam<Real>(
17 : "sat_lr",
18 22822 : 0.0,
19 : "sat_lr >= 0 & sat_lr < 1",
20 : "Liquid residual saturation. Must be between 0 and 1. Default is 0");
21 34233 : params.addRangeCheckedParam<Real>("pc_max",
22 22822 : 1.0e9,
23 : "pc_max > 0",
24 : "Maximum capillary pressure (Pa). Must be > 0. Default is 1e9");
25 22822 : params.addParam<bool>("log_extension",
26 22822 : 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 11411 : params.addClassDescription("Capillary pressure base class");
33 11411 : return params;
34 0 : }
35 :
36 6251 : PorousFlowCapillaryPressure::PorousFlowCapillaryPressure(const InputParameters & parameters)
37 : : DiscreteElementUserObject(parameters),
38 6251 : _sat_lr(getParam<Real>("sat_lr")),
39 6251 : _dseff_ds(1.0 / (1.0 - _sat_lr)),
40 12502 : _log_ext(getParam<bool>("log_extension")),
41 12502 : _pc_max(getParam<Real>("pc_max")),
42 6251 : _sat_ext(0.0),
43 6251 : _pc_ext(0.0),
44 6251 : _slope_ext(0.0),
45 6251 : _log10(std::log(10.0))
46 : {
47 6251 : }
48 :
49 : void
50 6020 : PorousFlowCapillaryPressure::initialSetup()
51 : {
52 : // If _log_ext = true, calculate the saturation where the the logarithmic
53 : // extension meets the raw capillary pressure curve.
54 6020 : if (_log_ext)
55 : {
56 3974 : _sat_ext = extensionSaturation();
57 3974 : _pc_ext = capillaryPressure(_sat_ext);
58 3974 : _slope_ext = (std::log10(_pc_ext) - std::log10(_pc_max)) / _sat_ext;
59 : }
60 6020 : }
61 :
62 : Real
63 4808132 : PorousFlowCapillaryPressure::capillaryPressure(Real saturation, unsigned qp) const
64 : {
65 4808132 : if (_log_ext && saturation < _sat_ext)
66 9782 : return capillaryPressureLogExt(saturation);
67 : else
68 4798350 : return capillaryPressureCurve(saturation, qp);
69 : }
70 :
71 : ADReal
72 1678597 : PorousFlowCapillaryPressure::capillaryPressure(const ADReal & saturation, unsigned qp) const
73 : {
74 1678597 : const Real Pc = capillaryPressure(saturation.value(), qp);
75 1678597 : const Real dPc_ds = dCapillaryPressure(saturation.value(), qp);
76 :
77 1678597 : ADReal result = Pc;
78 1678597 : result.derivatives() = saturation.derivatives() * dPc_ds;
79 :
80 1678597 : return result;
81 : }
82 :
83 : Real
84 4064334 : PorousFlowCapillaryPressure::dCapillaryPressure(Real saturation, unsigned qp) const
85 : {
86 4064334 : if (_log_ext && saturation < _sat_ext)
87 9038 : return dCapillaryPressureLogExt(saturation);
88 : else
89 4055296 : return dCapillaryPressureCurve(saturation, qp);
90 : }
91 :
92 : ADReal
93 34060 : PorousFlowCapillaryPressure::dCapillaryPressure(const ADReal & saturation, unsigned qp) const
94 : {
95 34060 : const Real dPc = dCapillaryPressure(saturation.value(), qp);
96 34060 : const Real d2Pc_ds2 = d2CapillaryPressure(saturation.value(), qp);
97 :
98 34060 : ADReal result = dPc;
99 34060 : result.derivatives() = saturation.derivatives() * d2Pc_ds2;
100 :
101 34060 : return result;
102 : }
103 :
104 : Real
105 1527761 : PorousFlowCapillaryPressure::d2CapillaryPressure(Real saturation, unsigned qp) const
106 : {
107 1527761 : if (_log_ext && saturation < _sat_ext)
108 7062 : return d2CapillaryPressureLogExt(saturation);
109 : else
110 1520699 : return d2CapillaryPressureCurve(saturation, qp);
111 : }
112 :
113 : Real
114 3313820 : PorousFlowCapillaryPressure::effectiveSaturationFromSaturation(Real saturation) const
115 : {
116 3313820 : return (saturation - _sat_lr) / (1.0 - _sat_lr);
117 : }
118 :
119 : Real
120 47852275 : PorousFlowCapillaryPressure::saturation(Real pc, unsigned qp) const
121 : {
122 47852275 : return effectiveSaturation(pc, qp) * (1.0 - _sat_lr) + _sat_lr;
123 : }
124 :
125 : ADReal
126 355173 : PorousFlowCapillaryPressure::saturation(const ADReal & pc, unsigned qp) const
127 : {
128 355173 : const Real s = effectiveSaturation(pc.value(), qp) * (1.0 - _sat_lr) + _sat_lr;
129 355173 : const Real ds_dpc = dEffectiveSaturation(pc.value(), qp) * (1.0 - _sat_lr);
130 :
131 355173 : ADReal result = s;
132 355173 : result.derivatives() = pc.derivatives() * ds_dpc;
133 :
134 355173 : return result;
135 : }
136 :
137 : Real
138 46651921 : PorousFlowCapillaryPressure::dSaturation(Real pc, unsigned qp) const
139 : {
140 46651921 : return dEffectiveSaturation(pc, qp) * (1.0 - _sat_lr);
141 : }
142 :
143 : ADReal
144 353642 : PorousFlowCapillaryPressure::dSaturation(const ADReal & pc, unsigned qp) const
145 : {
146 353642 : const Real ds = dEffectiveSaturation(pc.value(), qp) * (1.0 - _sat_lr);
147 353642 : const Real d2s_dpc2 = d2EffectiveSaturation(pc.value(), qp) * (1.0 - _sat_lr);
148 :
149 353642 : ADReal result = ds;
150 353642 : result.derivatives() = pc.derivatives() * d2s_dpc2;
151 :
152 353642 : return result;
153 : }
154 :
155 : Real
156 23137461 : PorousFlowCapillaryPressure::d2Saturation(Real pc, unsigned qp) const
157 : {
158 23137461 : return d2EffectiveSaturation(pc, qp) * (1.0 - _sat_lr);
159 : }
160 :
161 : Real
162 9782 : PorousFlowCapillaryPressure::capillaryPressureLogExt(Real saturation) const
163 : {
164 9782 : return _pc_ext * std::pow(10.0, _slope_ext * (saturation - _sat_ext));
165 : }
166 :
167 : Real
168 9038 : PorousFlowCapillaryPressure::dCapillaryPressureLogExt(Real saturation) const
169 : {
170 9038 : return _pc_ext * _slope_ext * _log10 * std::pow(10.0, _slope_ext * (saturation - _sat_ext));
171 : }
172 :
173 : Real
174 7062 : PorousFlowCapillaryPressure::d2CapillaryPressureLogExt(Real saturation) const
175 : {
176 7062 : return _pc_ext * _slope_ext * _slope_ext * _log10 * _log10 *
177 7062 : std::pow(10.0, _slope_ext * (saturation - _sat_ext));
178 : }
179 :
180 : Real
181 3974 : PorousFlowCapillaryPressure::extensionSaturation() const
182 : {
183 : // Initial guess for saturation where extension matches curve
184 3974 : 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 10368 : while (std::abs(interceptFunction(sat)) > nr_tol)
193 : {
194 6394 : sat = sat - interceptFunction(sat) / interceptFunctionDeriv(sat);
195 :
196 6394 : iter++;
197 6394 : if (iter > max_its)
198 : break;
199 : }
200 :
201 3974 : return sat;
202 : }
203 :
204 : Real
205 16762 : PorousFlowCapillaryPressure::interceptFunction(Real saturation) const
206 : {
207 16762 : Real pc = capillaryPressureCurve(saturation);
208 16762 : Real dpc = dCapillaryPressureCurve(saturation);
209 :
210 16762 : return std::log10(pc) - saturation * dpc / (_log10 * pc) - std::log10(_pc_max);
211 : }
212 :
213 : Real
214 6394 : PorousFlowCapillaryPressure::interceptFunctionDeriv(Real saturation) const
215 : {
216 6394 : Real pc = capillaryPressureCurve(saturation);
217 6394 : Real dpc = dCapillaryPressureCurve(saturation);
218 6394 : Real d2pc = d2CapillaryPressureCurve(saturation);
219 :
220 6394 : return saturation * (dpc * dpc / pc - d2pc) / (_log10 * pc);
221 : }
|
