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