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 "PorousFlowHystereticCapillaryPressure.h"
11 : #include "PorousFlowCapillaryPressure.h"
12 :
13 : registerMooseObject("PorousFlowApp", PorousFlowHystereticCapillaryPressure);
14 :
15 : InputParameters
16 10001 : PorousFlowHystereticCapillaryPressure::validParams()
17 : {
18 10001 : InputParameters params = PorousFlowVariableBase::validParams();
19 20002 : params.addRequiredRangeCheckedParam<Real>(
20 : "alpha_d",
21 : "alpha_d > 0",
22 : "van Genuchten alpha parameter for the primary drying curve. If using standard units, this "
23 : "is measured in Pa^-1. Suggested value is around 1E-5");
24 20002 : params.addRequiredRangeCheckedParam<Real>(
25 : "alpha_w",
26 : "alpha_w > 0",
27 : "van Genuchten alpha parameter for the primary wetting curve. If using standard units, this "
28 : "is measured in Pa^-1. Suggested value is around 1E-5");
29 20002 : params.addRequiredRangeCheckedParam<Real>("n_d",
30 : "n_d > 1",
31 : "van Genuchten n parameter for the primary drying "
32 : "curve. Dimensionless. Suggested value is around 2");
33 20002 : params.addRequiredRangeCheckedParam<Real>("n_w",
34 : "n_w > 1",
35 : "van Genuchten n parameter for the primary wetting "
36 : "curve. Dimensionless. Suggested value is around 2");
37 20002 : params.addRequiredRangeCheckedParam<Real>(
38 : "S_l_min",
39 : "S_l_min >= 0 & S_l_min < 1",
40 : "Minimum liquid saturation for which the van Genuchten expression is valid. If no lower "
41 : "extension is used then Pc = Pc_max for liquid saturation <= S_l_min");
42 30003 : params.addRangeCheckedParam<Real>(
43 : "S_lr",
44 20002 : 0.0,
45 : "S_lr >= 0 & S_lr < 1",
46 : "Liquid residual saturation where the liquid relative permeability is zero. This is used in "
47 : "the Land expression to find S_gr_del. Almost definitely you need to set S_lr equal to the "
48 : "quantity used for your relative-permeability curves. Almost definitely you should set S_lr "
49 : "> S_l_min");
50 20002 : params.addRequiredRangeCheckedParam<Real>(
51 : "S_gr_max",
52 : "S_gr_max >= 0",
53 : "Residual gas saturation. 1 - S_gr_max is the maximum saturation for which the "
54 : "van Genuchten expression is valid for the wetting curve. You must ensure S_gr_max < 1 - "
55 : "S_l_min. Often S_gr_max = -0.3136 * ln(porosity) - 0.1334 is used");
56 30003 : params.addRangeCheckedParam<Real>(
57 : "Pc_max",
58 20002 : std::numeric_limits<Real>::max(),
59 : "Pc_max > 0",
60 : "Value of capillary pressure at which the lower extension commences. The default value "
61 : "means capillary pressure uses the van Genuchten expression for S > S_l_min and is "
62 : "'infinity' for S <= S_l_min. This will result in poor convergence around S = S_l_min");
63 20002 : MooseEnum low_ext_enum("none quadratic exponential", "exponential");
64 20002 : params.addParam<MooseEnum>(
65 : "low_extension_type",
66 : low_ext_enum,
67 : "Type of extension to use for small liquid saturation values. The extensions modify the "
68 : "capillary pressure for all saturation values less than S(Pc_max). That is, if the "
69 : "van Genuchten "
70 : "expression would produce Pc > Pc_max, then the extension is used instead. NONE: Simply "
71 : "cut-off the capillary-pressure at Pc_max, so that Pc <= Pc_max for all S. QUADRATIC: Pc is "
72 : "a quadratic in S that is continuous and differentiable at S(Pc_max) and has zero derivative "
73 : "at S = 0 (hence, its value at S = 0 will be greater than Pc_max). EXPONENTIAL: Pc is an "
74 : "exponential in S that is continuous and differentiable at S(Pc_max) (hence, its value at S "
75 : "= 0 will be much greater than Pc_max");
76 30003 : params.addRangeCheckedParam<Real>(
77 : "high_ratio",
78 20002 : 0.9,
79 : "high_ratio > 0 & high_ratio < 1",
80 : "The extension to the wetting curves commences at high_ratio * (1 - S_gr_del), where 1 - "
81 : "S_gr_del is the value of the liquid saturation when Pc = 0 (on the wetting curve)");
82 20002 : MooseEnum high_ext_enum("none power", "power");
83 20002 : params.addParam<MooseEnum>(
84 : "high_extension_type",
85 : high_ext_enum,
86 : "Type of extension to use for the wetting curves when the liquid saturation is around 1. "
87 : "The extensions modify the wetting capillary pressure for all saturation values greater than "
88 : "high_ratio * (1 - S_gr_del), where 1 - S_gr_del is the value of liquid saturation when the "
89 : "van Genuchten expression gives Pc = 0. NONE: use the van Genuchten expression and when S > "
90 : "1 - S_gr_del, set Pc = 0. POWER: Pc is proportional to (1 - S)^power, where the "
91 : "coefficient of proportionality and the power are chosen so the resulting curve is "
92 : "continuous and differentiable");
93 10001 : params.addClassDescription(
94 : "This Material computes information that is required for the computation of hysteretic "
95 : "capillary pressures in single and multi phase situations");
96 10001 : return params;
97 10001 : }
98 :
99 7770 : PorousFlowHystereticCapillaryPressure::PorousFlowHystereticCapillaryPressure(
100 7770 : const InputParameters & parameters)
101 : : PorousFlowVariableBase(parameters),
102 7770 : _alpha_d(getParam<Real>("alpha_d")),
103 15540 : _alpha_w(getParam<Real>("alpha_w")),
104 15540 : _n_d(getParam<Real>("n_d")),
105 15540 : _n_w(getParam<Real>("n_w")),
106 15540 : _s_l_min(getParam<Real>("S_l_min")),
107 15540 : _s_lr(getParam<Real>("S_lr")),
108 15540 : _s_gr_max(getParam<Real>("S_gr_max")),
109 15540 : _pc_max(getParam<Real>("Pc_max")),
110 15540 : _high_ratio(getParam<Real>("high_ratio")),
111 7770 : _low_ext_type(
112 15540 : getParam<MooseEnum>("low_extension_type")
113 : .getEnum<PorousFlowVanGenuchten::LowCapillaryPressureExtension::ExtensionStrategy>()),
114 7770 : _s_low_d(PorousFlowVanGenuchten::saturationHys(_pc_max, _s_l_min, 0.0, _alpha_d, _n_d)),
115 7770 : _dpc_low_d(
116 7770 : PorousFlowVanGenuchten::dcapillaryPressureHys(_s_low_d, _s_l_min, 0.0, _alpha_d, _n_d)),
117 7770 : _low_ext_d(_low_ext_type, _s_low_d, _pc_max, _dpc_low_d),
118 7770 : _s_low_w(PorousFlowVanGenuchten::saturationHys(_pc_max, _s_l_min, _s_gr_max, _alpha_w, _n_w)),
119 15540 : _dpc_low_w(PorousFlowVanGenuchten::dcapillaryPressureHys(
120 7770 : _s_low_w, _s_l_min, _s_gr_max, _alpha_w, _n_w)),
121 7770 : _low_ext_w(_low_ext_type, _s_low_w, _pc_max, _dpc_low_w),
122 7770 : _high_ext_type(
123 7770 : getParam<MooseEnum>("high_extension_type")
124 : .getEnum<PorousFlowVanGenuchten::HighCapillaryPressureExtension::ExtensionStrategy>()),
125 7770 : _s_high(_high_ratio * (1 - _s_gr_max)),
126 7770 : _pc_high(
127 7770 : PorousFlowVanGenuchten::capillaryPressureHys(_s_high, _s_l_min, _s_gr_max, _alpha_w, _n_w)),
128 15540 : _dpc_high(PorousFlowVanGenuchten::dcapillaryPressureHys(
129 7770 : _s_high, _s_l_min, _s_gr_max, _alpha_w, _n_w)),
130 7770 : _high_ext(_high_ext_type, _s_high, _pc_high, _dpc_high),
131 7770 : _hys_order(_nodal_material ? getMaterialProperty<unsigned>("PorousFlow_hysteresis_order_nodal")
132 22342 : : getMaterialProperty<unsigned>("PorousFlow_hysteresis_order_qp")),
133 15540 : _hys_order_old(_nodal_material
134 7770 : ? getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_nodal")
135 22342 : : getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_qp")),
136 7770 : _hys_sat_tps(
137 7770 : _nodal_material
138 7770 : ? getMaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
139 : "PorousFlow_hysteresis_saturation_tps_nodal")
140 22342 : : getMaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
141 : "PorousFlow_hysteresis_saturation_tps_qp")),
142 15540 : _pc_older(_nodal_material
143 7770 : ? getMaterialPropertyOlder<Real>("PorousFlow_hysteretic_capillary_pressure_nodal")
144 22342 : : getMaterialPropertyOlder<Real>("PorousFlow_hysteretic_capillary_pressure_qp")),
145 15540 : _pc_tps(_nodal_material
146 7770 : ? declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
147 : "PorousFlow_hysteresis_Pc_tps_nodal")
148 15056 : : declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
149 : "PorousFlow_hysteresis_Pc_tps_qp")),
150 15540 : _s_d_tps(_nodal_material
151 7770 : ? declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
152 : "PorousFlow_hysteresis_s_d_tps_nodal")
153 15056 : : declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
154 : "PorousFlow_hysteresis_s_d_tps_qp")),
155 15540 : _s_gr_tps(_nodal_material
156 7770 : ? declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
157 : "PorousFlow_hysteresis_s_gr_tps_nodal")
158 15056 : : declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
159 : "PorousFlow_hysteresis_s_gr_tps_qp")),
160 7770 : _w_low_ext_tps(
161 7770 : _nodal_material
162 7770 : ? declareProperty<std::array<PorousFlowVanGenuchten::LowCapillaryPressureExtension,
163 484 : PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
164 : "PorousFlow_hysteresis_w_low_ext_tps_nodal")
165 : : declareProperty<std::array<PorousFlowVanGenuchten::LowCapillaryPressureExtension,
166 15056 : PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
167 : "PorousFlow_hysteresis_w_low_ext_tps_qp")),
168 7770 : _w_high_ext_tps(
169 7770 : _nodal_material
170 7770 : ? declareProperty<std::array<PorousFlowVanGenuchten::HighCapillaryPressureExtension,
171 484 : PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
172 : "PorousFlow_hysteresis_w_high_ext_tps_nodal")
173 : : declareProperty<std::array<PorousFlowVanGenuchten::HighCapillaryPressureExtension,
174 15056 : PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
175 : "PorousFlow_hysteresis_w_high_ext_tps_qp")),
176 15540 : _s_w_tps(_nodal_material
177 7770 : ? declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
178 : "PorousFlow_hysteresis_s_w_tps_nodal")
179 15056 : : declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
180 7770 : "PorousFlow_hysteresis_s_w_tps_qp"))
181 : {
182 7770 : if (_s_gr_max >= 1 - _s_l_min)
183 2 : paramError("S_gr_max", "Must be less than 1 - S_l_min");
184 7768 : if (_s_high < _s_low_w)
185 2 : paramError("high_ratio",
186 : "Should be chosen sufficiently close to 1 so that the high extension does not "
187 : "interfere with the low extension. Instead, you may have chosen Pc_max too low");
188 7766 : if (_s_lr <= _s_l_min)
189 2 : mooseWarning("S_lr should usually be greater than S_l_min");
190 7764 : }
191 :
192 : void
193 344096 : PorousFlowHystereticCapillaryPressure::initQpStatefulProperties()
194 : {
195 344096 : PorousFlowVariableBase::initQpStatefulProperties();
196 : Real tp_sat;
197 : Real tp_pc;
198 344096 : if (_hys_order[_qp] >= 1)
199 : {
200 254988 : tp_sat = _hys_sat_tps[_qp].at(0);
201 509976 : tp_pc = PorousFlowVanGenuchten::capillaryPressureHys(
202 254988 : tp_sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
203 254988 : computeTurningPointInfo(0, tp_sat, tp_pc);
204 : }
205 344096 : if (_hys_order[_qp] >= 2)
206 : {
207 128940 : tp_sat = _hys_sat_tps[_qp].at(1);
208 128940 : tp_pc = firstOrderWettingPc(tp_sat);
209 128940 : computeTurningPointInfo(1, tp_sat, tp_pc);
210 : }
211 344096 : if (_hys_order[_qp] >= 3)
212 : {
213 48300 : tp_sat = _hys_sat_tps[_qp].at(2);
214 48300 : tp_pc = secondOrderDryingPc(tp_sat);
215 48300 : computeTurningPointInfo(2, tp_sat, tp_pc);
216 : }
217 344096 : }
218 :
219 : void
220 824538 : PorousFlowHystereticCapillaryPressure::computeQpProperties()
221 : {
222 : // size stuff correctly and prepare the derivative matrices with zeroes
223 824538 : PorousFlowVariableBase::computeQpProperties();
224 :
225 824538 : if (_hys_order[_qp] != _hys_order_old[_qp] && _hys_order[_qp] > 0)
226 : {
227 14478 : const unsigned tp_num = _hys_order[_qp] - 1;
228 14478 : computeTurningPointInfo(tp_num, _hys_sat_tps[_qp].at(tp_num), _pc_older[_qp]);
229 : }
230 824538 : }
231 :
232 : Real
233 446706 : PorousFlowHystereticCapillaryPressure::landSat(Real slDel) const
234 : {
235 446706 : const Real a = 1.0 / _s_gr_max - 1.0 / (1.0 - _s_lr);
236 446706 : return (1.0 - slDel) / (1.0 + a * (1.0 - slDel));
237 : }
238 :
239 : void
240 446706 : PorousFlowHystereticCapillaryPressure::computeTurningPointInfo(unsigned tp_num,
241 : Real tp_sat,
242 : Real tp_pc)
243 : {
244 446706 : _pc_tps[_qp].at(tp_num) = tp_pc;
245 :
246 : // Quantities on the drying curve:
247 : // pc on the drying curve at the turning point saturation
248 446706 : const Real pc_d_tps = PorousFlowVanGenuchten::capillaryPressureHys(
249 446706 : tp_sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
250 : // saturation on the drying curve, at tp_pc
251 446706 : _s_d_tps[_qp].at(tp_num) =
252 446706 : PorousFlowVanGenuchten::saturationHys(tp_pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
253 :
254 : // Quantities relevant to the wetting curve defined by the Land expression.
255 : // s_gr_tps is the Land expression as a function of the turning point saturation
256 446706 : _s_gr_tps[_qp].at(tp_num) = landSat(tp_sat);
257 : // the low extension of the wetting curve defined by _s_gr_tps
258 446706 : const Real s_w_low_ext = PorousFlowVanGenuchten::saturationHys(
259 446706 : _pc_max, _s_l_min, _s_gr_tps[_qp].at(tp_num), _alpha_w, _n_w);
260 446706 : const Real dpc_w_low_ext = PorousFlowVanGenuchten::dcapillaryPressureHys(
261 446706 : s_w_low_ext, _s_l_min, _s_gr_tps[_qp].at(tp_num), _alpha_w, _n_w);
262 446706 : _w_low_ext_tps[_qp].at(tp_num) = PorousFlowVanGenuchten::LowCapillaryPressureExtension(
263 446706 : _low_ext_type, s_w_low_ext, _pc_max, dpc_w_low_ext);
264 : // the high extension of the wetting curve defined by _s_gr_tps
265 : const Real s_w_high_ext =
266 446706 : (_high_ext_type == PorousFlowVanGenuchten::HighCapillaryPressureExtension::NONE)
267 446706 : ? 1.0 - _s_gr_tps[_qp].at(tp_num)
268 268206 : : _high_ratio *
269 268206 : (1.0 -
270 268206 : _s_gr_tps[_qp].at(
271 : tp_num)); // if NONE then use the vanGenuchten all the way to 1 - _s_gr_tps
272 : const Real pc_w_high_ext =
273 446706 : PorousFlowVanGenuchten::capillaryPressureHys(s_w_high_ext,
274 446706 : _s_l_min,
275 446706 : _s_gr_tps[_qp].at(tp_num),
276 446706 : _alpha_w,
277 446706 : _n_w,
278 : _w_low_ext_tps[_qp].at(tp_num));
279 : const Real dpc_w_high_ext =
280 446706 : PorousFlowVanGenuchten::dcapillaryPressureHys(s_w_high_ext,
281 446706 : _s_l_min,
282 446706 : _s_gr_tps[_qp].at(tp_num),
283 446706 : _alpha_w,
284 446706 : _n_w,
285 446706 : _w_low_ext_tps[_qp].at(tp_num));
286 446706 : _w_high_ext_tps[_qp].at(tp_num) = PorousFlowVanGenuchten::HighCapillaryPressureExtension(
287 : _high_ext_type, s_w_high_ext, pc_w_high_ext, dpc_w_high_ext);
288 :
289 : // saturation on the wetting curve defined by _s_gr_tps, at pc = pc_d_tps
290 893412 : _s_w_tps[_qp].at(tp_num) = PorousFlowVanGenuchten::saturationHys(pc_d_tps,
291 446706 : _s_l_min,
292 446706 : _s_gr_tps[_qp].at(tp_num),
293 446706 : _alpha_w,
294 446706 : _n_w,
295 446706 : _w_low_ext_tps[_qp].at(tp_num),
296 : _w_high_ext_tps[_qp].at(tp_num));
297 446706 : }
298 :
299 : Real
300 1966000 : PorousFlowHystereticCapillaryPressure::capillaryPressureQp(Real sat) const
301 : {
302 : Real pc = 0.0;
303 1966000 : if (_hys_order[_qp] == 0) // on primary drying curve
304 464596 : pc = PorousFlowVanGenuchten::capillaryPressureHys(
305 464596 : sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
306 1501404 : else if (_hys_order[_qp] == 1) // first-order wetting
307 746348 : pc = firstOrderWettingPc(sat);
308 755056 : else if (_hys_order[_qp] == 2) // second-order drying
309 473016 : pc = secondOrderDryingPc(sat);
310 : else // third order drying and wetting
311 : {
312 282040 : const Real tp1 = _hys_sat_tps[_qp].at(1);
313 282040 : const Real tp2 = _hys_sat_tps[_qp].at(2);
314 282040 : const Real pc1 = firstOrderWettingPc(sat);
315 282040 : const Real pc2 = secondOrderDryingPc(sat);
316 : // handle cases that occur just at the transition from 3rd to 2nd order, or 3rd to 1st order
317 282040 : if (sat < tp2)
318 : pc = pc2;
319 282040 : else if (sat > tp1)
320 : pc = pc1;
321 281896 : else if (pc1 >= pc2)
322 : pc = pc1;
323 281896 : else if (pc1 <= 0.0 || pc2 <= 0.0)
324 : pc = 0.0;
325 : else
326 234682 : pc = std::exp(std::log(pc1) + (sat - tp1) * (std::log(pc2) - std::log(pc1)) / (tp2 - tp1));
327 : }
328 1966000 : return pc;
329 : }
330 :
331 : Real
332 389808 : PorousFlowHystereticCapillaryPressure::dcapillaryPressureQp(Real sat) const
333 : {
334 : Real dpc = 0.0;
335 389808 : if (_hys_order[_qp] == 0) // on primary drying curve
336 125652 : dpc = PorousFlowVanGenuchten::dcapillaryPressureHys(
337 125652 : sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
338 264156 : else if (_hys_order[_qp] == 1) // first-order wetting
339 119040 : dpc = dfirstOrderWettingPc(sat);
340 145116 : else if (_hys_order[_qp] == 2) // second-order drying
341 90216 : dpc = dsecondOrderDryingPc(sat);
342 : else // third order drying and wetting
343 : {
344 54900 : const Real tp1 = _hys_sat_tps[_qp].at(1);
345 54900 : const Real tp2 = _hys_sat_tps[_qp].at(2);
346 54900 : const Real pc1 = firstOrderWettingPc(sat);
347 54900 : const Real pc2 = secondOrderDryingPc(sat);
348 : // handle cases that occur just at the transition from 3rd to 2nd order, or 3rd to 1st order
349 54900 : if (sat < tp2)
350 0 : dpc = dsecondOrderDryingPc(sat);
351 54900 : else if (sat > tp1)
352 0 : dpc = dfirstOrderWettingPc(sat);
353 54900 : else if (pc1 >= pc2)
354 0 : dpc = dfirstOrderWettingPc(sat);
355 54900 : else if (pc1 <= 0.0 || pc2 <= 0.0)
356 : dpc = 0.0;
357 : else
358 : {
359 : const Real pc =
360 47031 : std::exp(std::log(pc1) + (sat - tp1) * (std::log(pc2) - std::log(pc1)) / (tp2 - tp1));
361 47031 : const Real dpc1 = dfirstOrderWettingPc(sat);
362 47031 : const Real dpc2 = dsecondOrderDryingPc(sat);
363 47031 : dpc = pc * (dpc1 / pc1 + (std::log(pc2) - std::log(pc1)) / (tp2 - tp1) +
364 47031 : (sat - tp1) * (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1));
365 : }
366 : }
367 389808 : return dpc;
368 : }
369 :
370 : Real
371 134514 : PorousFlowHystereticCapillaryPressure::d2capillaryPressureQp(Real sat) const
372 : {
373 : Real d2pc = 0.0;
374 134514 : if (_hys_order[_qp] == 0) // on primary drying curve
375 44526 : d2pc = PorousFlowVanGenuchten::d2capillaryPressureHys(
376 44526 : sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
377 89988 : else if (_hys_order[_qp] == 1) // first-order wetting
378 41220 : d2pc = d2firstOrderWettingPc(sat);
379 48768 : else if (_hys_order[_qp] == 2) // second-order drying
380 30468 : d2pc = d2secondOrderDryingPc(sat);
381 : else // third order drying and wetting
382 : {
383 18300 : const Real tp1 = _hys_sat_tps[_qp].at(1);
384 18300 : const Real tp2 = _hys_sat_tps[_qp].at(2);
385 18300 : const Real pc1 = firstOrderWettingPc(sat);
386 18300 : const Real pc2 = secondOrderDryingPc(sat);
387 : // handle cases that occur just at the transition from 3rd to 2nd order, or 3rd to 1st order
388 18300 : if (sat < tp2)
389 0 : d2pc = d2secondOrderDryingPc(sat);
390 18300 : else if (sat > tp1)
391 0 : d2pc = d2firstOrderWettingPc(sat);
392 18300 : else if (pc1 >= pc2)
393 0 : d2pc = d2firstOrderWettingPc(sat);
394 18300 : else if (pc1 <= 0.0 || pc2 <= 0.0)
395 : d2pc = 0.0;
396 : else
397 : {
398 : const Real pc =
399 15677 : std::exp(std::log(pc1) + (sat - tp1) * (std::log(pc2) - std::log(pc1)) / (tp2 - tp1));
400 15677 : const Real dpc1 = dfirstOrderWettingPc(sat);
401 15677 : const Real dpc2 = dsecondOrderDryingPc(sat);
402 15677 : const Real dpc = pc * (dpc1 / pc1 + (std::log(pc2) - std::log(pc1)) / (tp2 - tp1) +
403 15677 : (sat - tp1) * (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1));
404 15677 : const Real d2pc1 = d2firstOrderWettingPc(sat);
405 15677 : const Real d2pc2 = d2secondOrderDryingPc(sat);
406 15677 : d2pc =
407 15677 : dpc * (dpc1 / pc1 + (std::log(pc2) - std::log(pc1)) / (tp2 - tp1) +
408 : (sat - tp1) * (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1)) +
409 15677 : pc *
410 15677 : (d2pc1 / pc1 - dpc1 * dpc1 / pc1 / pc1 + (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1) +
411 15677 : (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1) +
412 15677 : (sat - tp1) *
413 15677 : (d2pc2 / pc2 - dpc2 * dpc2 / pc2 / pc2 - d2pc1 / pc1 + dpc1 * dpc1 / pc1 / pc1) /
414 : (tp2 - tp1));
415 : }
416 : }
417 134514 : return d2pc;
418 : }
419 :
420 : Real
421 1230528 : PorousFlowHystereticCapillaryPressure::firstOrderWettingPc(Real sat) const
422 : {
423 : // Simplest version is to just use the wetting curve defined by _s_gr_tps[0], but want
424 : // continuity at sat = _hys_sat_tps[0] (the turning point), so use the following process. The
425 : // wetting curve is defined for S <= max_s, where
426 1230528 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
427 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
428 1230528 : ? 1.0
429 503900 : : 1.0 - _s_gr_tps[_qp].at(0);
430 : // define an interpolation: s_to_use smoothly transitions from _s_w_tps[0] (the value of liquid
431 : // saturation on the wetting curve defined by _s_gr_tps[0]) when sat = _hys_sat_tps[0], to max_s
432 : // when sat = max_s
433 1230528 : const Real sat_to_use = _s_w_tps[_qp].at(0) + (max_s - _s_w_tps[_qp].at(0)) *
434 1230528 : (sat - _hys_sat_tps[_qp].at(0)) /
435 1230528 : (max_s - _hys_sat_tps[_qp].at(0));
436 1230528 : return PorousFlowVanGenuchten::capillaryPressureHys(sat_to_use,
437 1230528 : _s_l_min,
438 1230528 : _s_gr_tps[_qp].at(0),
439 1230528 : _alpha_w,
440 1230528 : _n_w,
441 1230528 : _w_low_ext_tps[_qp].at(0),
442 1230528 : _w_high_ext_tps[_qp].at(0));
443 : }
444 :
445 : Real
446 181748 : PorousFlowHystereticCapillaryPressure::dfirstOrderWettingPc(Real sat) const
447 : {
448 181748 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
449 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
450 181748 : ? 1.0
451 50508 : : 1.0 - _s_gr_tps[_qp].at(0);
452 181748 : const Real sat_to_use = _s_w_tps[_qp].at(0) + (max_s - _s_w_tps[_qp].at(0)) *
453 181748 : (sat - _hys_sat_tps[_qp].at(0)) /
454 181748 : (max_s - _hys_sat_tps[_qp].at(0));
455 181748 : const Real dsat_to_use = (max_s - _s_w_tps[_qp].at(0)) / (max_s - _hys_sat_tps[_qp].at(0));
456 181748 : return PorousFlowVanGenuchten::dcapillaryPressureHys(sat_to_use,
457 181748 : _s_l_min,
458 181748 : _s_gr_tps[_qp].at(0),
459 181748 : _alpha_w,
460 181748 : _n_w,
461 181748 : _w_low_ext_tps[_qp].at(0),
462 : _w_high_ext_tps[_qp].at(0)) *
463 181748 : dsat_to_use;
464 : }
465 :
466 : Real
467 56897 : PorousFlowHystereticCapillaryPressure::d2firstOrderWettingPc(Real sat) const
468 : {
469 56897 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
470 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
471 56897 : ? 1.0
472 15677 : : 1.0 - _s_gr_tps[_qp].at(0);
473 56897 : const Real sat_to_use = _s_w_tps[_qp].at(0) + (max_s - _s_w_tps[_qp].at(0)) *
474 56897 : (sat - _hys_sat_tps[_qp].at(0)) /
475 56897 : (max_s - _hys_sat_tps[_qp].at(0));
476 56897 : const Real dsat_to_use = (max_s - _s_w_tps[_qp].at(0)) / (max_s - _hys_sat_tps[_qp].at(0));
477 56897 : return PorousFlowVanGenuchten::d2capillaryPressureHys(sat_to_use,
478 56897 : _s_l_min,
479 56897 : _s_gr_tps[_qp].at(0),
480 56897 : _alpha_w,
481 56897 : _n_w,
482 56897 : _w_low_ext_tps[_qp].at(0),
483 56897 : _w_high_ext_tps[_qp].at(0)) *
484 56897 : dsat_to_use * dsat_to_use;
485 : }
486 :
487 : Real
488 876556 : PorousFlowHystereticCapillaryPressure::secondOrderDryingPc(Real sat) const
489 : {
490 : // Simplest version is to just use the primary drying curve, but want
491 : // continuity at sat = _hys_sat_tps[0] (the dry-to-wet turning point) and sat = _hys_sat_tps[1]
492 : // (the wet-to-dry turning point), so use the following process.
493 876556 : const Real tp0 = _hys_sat_tps[_qp].at(0);
494 876556 : const Real tp1 = _hys_sat_tps[_qp].at(1);
495 876556 : const Real s1 = _s_d_tps[_qp].at(1);
496 : const Real sat_to_use =
497 876556 : (sat >= tp0) ? tp0 + (sat - tp0) * (s1 - tp0) / (tp1 - tp0)
498 : : sat; // final case can occur just at transition from 2nd to 0th order
499 1753112 : return PorousFlowVanGenuchten::capillaryPressureHys(
500 876556 : sat_to_use, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
501 : }
502 :
503 : Real
504 152924 : PorousFlowHystereticCapillaryPressure::dsecondOrderDryingPc(Real sat) const
505 : {
506 152924 : const Real tp0 = _hys_sat_tps[_qp].at(0);
507 152924 : const Real tp1 = _hys_sat_tps[_qp].at(1);
508 152924 : const Real s1 = _s_d_tps[_qp].at(1);
509 152924 : const Real sat_to_use = (sat >= tp0) ? tp0 + (sat - tp0) * (s1 - tp0) / (tp1 - tp0) : sat;
510 152924 : const Real dsat_to_use = (sat >= tp0) ? (s1 - tp0) / (tp1 - tp0) : 1.0;
511 305848 : return PorousFlowVanGenuchten::dcapillaryPressureHys(
512 152924 : sat_to_use, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d) *
513 152924 : dsat_to_use;
514 : }
515 :
516 : Real
517 46145 : PorousFlowHystereticCapillaryPressure::d2secondOrderDryingPc(Real sat) const
518 : {
519 46145 : const Real tp0 = _hys_sat_tps[_qp].at(0);
520 46145 : const Real tp1 = _hys_sat_tps[_qp].at(1);
521 46145 : const Real s1 = _s_d_tps[_qp].at(1);
522 46145 : const Real sat_to_use = (sat >= tp0) ? tp0 + (sat - tp0) * (s1 - tp0) / (tp1 - tp0) : sat;
523 46145 : const Real dsat_to_use = (sat >= tp0) ? (s1 - tp0) / (tp1 - tp0) : 1.0;
524 92290 : return PorousFlowVanGenuchten::d2capillaryPressureHys(
525 46145 : sat_to_use, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d) *
526 46145 : dsat_to_use * dsat_to_use;
527 : }
528 :
529 : Real
530 358598 : PorousFlowHystereticCapillaryPressure::firstOrderWettingSat(Real pc) const
531 : {
532 : // this is inverse of firstOrderWettingPc: see that method for comments
533 358598 : const Real sat_to_use = PorousFlowVanGenuchten::saturationHys(pc,
534 358598 : _s_l_min,
535 358598 : _s_gr_tps[_qp].at(0),
536 358598 : _alpha_w,
537 358598 : _n_w,
538 358598 : _w_low_ext_tps[_qp].at(0),
539 358598 : _w_high_ext_tps[_qp].at(0));
540 358598 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
541 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
542 358598 : ? 1.0
543 109800 : : 1.0 - _s_gr_tps[_qp].at(0);
544 358598 : if (sat_to_use > max_s) // this occurs when using no high extension and pc = 0
545 : return max_s;
546 327427 : return (sat_to_use - _s_w_tps[_qp].at(0)) * (max_s - _hys_sat_tps[_qp].at(0)) /
547 327427 : (max_s - _s_w_tps[_qp].at(0)) +
548 327427 : _hys_sat_tps[_qp].at(0);
549 : }
550 :
551 : Real
552 247890 : PorousFlowHystereticCapillaryPressure::dfirstOrderWettingSat(Real pc) const
553 : {
554 247890 : const Real dsat_to_use = PorousFlowVanGenuchten::dsaturationHys(pc,
555 247890 : _s_l_min,
556 247890 : _s_gr_tps[_qp].at(0),
557 247890 : _alpha_w,
558 247890 : _n_w,
559 247890 : _w_low_ext_tps[_qp].at(0),
560 247890 : _w_high_ext_tps[_qp].at(0));
561 247890 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
562 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
563 247890 : ? 1.0
564 61000 : : 1.0 - _s_gr_tps[_qp].at(0);
565 247890 : if (pc <= 0.0 && _w_high_ext_tps[_qp].at(0).strategy ==
566 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::NONE)
567 : return 0.0;
568 232213 : return dsat_to_use * (max_s - _hys_sat_tps[_qp].at(0)) / (max_s - _s_w_tps[_qp].at(0));
569 : }
570 :
571 : Real
572 84510 : PorousFlowHystereticCapillaryPressure::d2firstOrderWettingSat(Real pc) const
573 : {
574 84510 : const Real d2sat_to_use = PorousFlowVanGenuchten::d2saturationHys(pc,
575 84510 : _s_l_min,
576 84510 : _s_gr_tps[_qp].at(0),
577 84510 : _alpha_w,
578 84510 : _n_w,
579 84510 : _w_low_ext_tps[_qp].at(0),
580 84510 : _w_high_ext_tps[_qp].at(0));
581 84510 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
582 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
583 84510 : ? 1.0
584 18300 : : 1.0 - _s_gr_tps[_qp].at(0);
585 84510 : if (pc <= 0.0 && _w_high_ext_tps[_qp].at(0).strategy ==
586 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::NONE)
587 : return 0.0;
588 77983 : return d2sat_to_use * (max_s - _hys_sat_tps[_qp].at(0)) / (max_s - _s_w_tps[_qp].at(0));
589 : }
590 :
591 : Real
592 285716 : PorousFlowHystereticCapillaryPressure::secondOrderDryingSat(Real pc) const
593 : {
594 : // This is the inverse of secondOrderDryingPc: see that method for comments
595 : const Real sat_to_use =
596 285716 : PorousFlowVanGenuchten::saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
597 285716 : const Real tp0 = _hys_sat_tps[_qp].at(0);
598 285716 : const Real tp1 = _hys_sat_tps[_qp].at(1);
599 285716 : const Real s1 = _s_d_tps[_qp].at(1);
600 285716 : return (sat_to_use >= tp0) ? (sat_to_use - tp0) * (tp1 - tp0) / (s1 - tp0) + tp0 : sat_to_use;
601 : }
602 :
603 : Real
604 182352 : PorousFlowHystereticCapillaryPressure::dsecondOrderDryingSat(Real pc) const
605 : {
606 : const Real sat_to_use =
607 182352 : PorousFlowVanGenuchten::saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
608 : const Real dsat_to_use =
609 182352 : PorousFlowVanGenuchten::dsaturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
610 182352 : const Real tp0 = _hys_sat_tps[_qp].at(0);
611 182352 : const Real tp1 = _hys_sat_tps[_qp].at(1);
612 182352 : const Real s1 = _s_d_tps[_qp].at(1);
613 182352 : return (sat_to_use >= tp0) ? dsat_to_use * (tp1 - tp0) / (s1 - tp0) : dsat_to_use;
614 : }
615 :
616 : Real
617 56022 : PorousFlowHystereticCapillaryPressure::d2secondOrderDryingSat(Real pc) const
618 : {
619 : const Real sat_to_use =
620 56022 : PorousFlowVanGenuchten::saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
621 : const Real d2sat_to_use =
622 56022 : PorousFlowVanGenuchten::d2saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
623 56022 : const Real tp0 = _hys_sat_tps[_qp].at(0);
624 56022 : const Real tp1 = _hys_sat_tps[_qp].at(1);
625 56022 : const Real s1 = _s_d_tps[_qp].at(1);
626 56022 : return (sat_to_use >= tp0) ? d2sat_to_use * (tp1 - tp0) / (s1 - tp0) : d2sat_to_use;
627 : }
628 :
629 : Real
630 612228 : PorousFlowHystereticCapillaryPressure::liquidSaturationQp(Real pc) const
631 : {
632 : Real sat = 0.0;
633 612228 : if (_hys_order[_qp] == 0) // on primary drying curve
634 209350 : sat = PorousFlowVanGenuchten::saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
635 402878 : else if (_hys_order[_qp] == 1) // first-order wetting
636 200838 : sat = firstOrderWettingSat(pc);
637 202040 : else if (_hys_order[_qp] == 2) // second-order drying
638 127956 : sat = secondOrderDryingSat(pc);
639 : else // third order drying and wetting
640 : {
641 : // NOTE: this is not the exact inverse of the third-order capillary-pressure formula, but only
642 : // the approximate inverse. In any simulation, only liquidSaturationQp or capillaryPressureQp
643 : // are used (not both) so having a slightly different formulation for these two functions is OK
644 : // Rationale: when pc is close to pc_tp1 then use the first-order wetting curve; when pc is
645 : // close to pc_tp2 then use the second-order drying curve
646 74084 : const Real pc_tp1 = _pc_tps[_qp].at(1); // pc on first-order wetting at TP_1
647 74084 : const Real pc_tp2 = _pc_tps[_qp].at(2); // pc on second-order drying at TP_2
648 74084 : const Real sat1 = firstOrderWettingSat(pc);
649 74084 : const Real sat2 = secondOrderDryingSat(pc);
650 : // handle cases that occur just at the transition from 3rd to 2nd order, or 3rd to 1st order
651 74084 : if (pc_tp1 <= 0.0 || pc <= 0.0 ||
652 : pc_tp2 <= 0.0) // only the first condition is strictly necessary as cannot get pc==0 or
653 : // pc_tp2=0 without pc_tp1=0 in reality. The other conditions are added in
654 : // case of numerical strangenesses
655 : sat = sat2;
656 49612 : else if (pc > pc_tp2)
657 : sat = sat2;
658 49612 : else if (pc < pc_tp1)
659 : sat = sat1;
660 : else
661 49072 : sat = sat1 + (std::log(pc) - std::log(pc_tp1)) * (sat2 - sat1) / std::log(pc_tp2 / pc_tp1);
662 : }
663 612228 : return sat;
664 : }
665 :
666 : Real
667 520986 : PorousFlowHystereticCapillaryPressure::dliquidSaturationQp(Real pc) const
668 : {
669 : Real dsat = 0.0;
670 520986 : if (_hys_order[_qp] == 0) // on primary drying curve
671 196212 : dsat = PorousFlowVanGenuchten::dsaturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
672 324774 : else if (_hys_order[_qp] == 1) // first-order wetting
673 164214 : dsat = dfirstOrderWettingSat(pc);
674 160560 : else if (_hys_order[_qp] == 2) // second-order drying
675 98676 : dsat = dsecondOrderDryingSat(pc);
676 : else // third order drying and wetting
677 : {
678 61884 : const Real pc_tp1 = _pc_tps[_qp].at(1); // pc on first-order wetting at TP_1
679 61884 : const Real pc_tp2 = _pc_tps[_qp].at(2); // pc on second-order drying at TP_2
680 61884 : const Real sat1 = firstOrderWettingSat(pc);
681 61884 : const Real sat2 = secondOrderDryingSat(pc);
682 61884 : const Real dsat1 = dfirstOrderWettingSat(pc);
683 61884 : const Real dsat2 = dsecondOrderDryingSat(pc);
684 61884 : if (pc_tp1 <= 0.0 || pc <= 0.0 || pc_tp2 <= 0.0)
685 : dsat = dsat2;
686 43512 : else if (pc > pc_tp2)
687 : dsat = dsat2;
688 43512 : else if (pc < pc_tp1)
689 : dsat = dsat1;
690 : else
691 42972 : dsat = dsat1 + 1.0 / pc * (sat2 - sat1) / std::log(pc_tp2 / pc_tp1) +
692 42972 : (std::log(pc) - std::log(pc_tp1)) * (dsat2 - dsat1) / std::log(pc_tp2 / pc_tp1);
693 : }
694 520986 : return dsat;
695 : }
696 :
697 : Real
698 208638 : PorousFlowHystereticCapillaryPressure::d2liquidSaturationQp(Real pc) const
699 : {
700 : Real d2sat = 0.0;
701 208638 : if (_hys_order[_qp] == 0) // on primary drying curve
702 89898 : d2sat = PorousFlowVanGenuchten::d2saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
703 118740 : else if (_hys_order[_qp] == 1) // first-order wetting
704 62718 : d2sat = d2firstOrderWettingSat(pc);
705 56022 : else if (_hys_order[_qp] == 2) // second-order drying
706 34230 : d2sat = d2secondOrderDryingSat(pc);
707 : else // third order drying and wetting
708 : {
709 21792 : const Real pc_tp1 = _pc_tps[_qp].at(1); // pc on first-order wetting at TP_1
710 21792 : const Real pc_tp2 = _pc_tps[_qp].at(2); // pc on second-order drying at TP_2
711 21792 : const Real sat1 = firstOrderWettingSat(pc);
712 21792 : const Real sat2 = secondOrderDryingSat(pc);
713 21792 : const Real dsat1 = dfirstOrderWettingSat(pc);
714 21792 : const Real dsat2 = dsecondOrderDryingSat(pc);
715 21792 : const Real d2sat1 = d2firstOrderWettingSat(pc);
716 21792 : const Real d2sat2 = d2secondOrderDryingSat(pc);
717 21792 : if (pc_tp1 <= 0.0 || pc <= 0.0 || pc_tp2 <= 0.0)
718 : d2sat = d2sat2;
719 15656 : else if (pc > pc_tp2)
720 : d2sat = d2sat2;
721 15656 : else if (pc < pc_tp1)
722 : d2sat = d2sat1;
723 : else
724 15386 : d2sat = d2sat1 - 1.0 / pc / pc * (sat2 - sat1) / std::log(pc_tp2 / pc_tp1) +
725 15386 : 2.0 / pc * (dsat2 - dsat1) / std::log(pc_tp2 / pc_tp1) +
726 15386 : (std::log(pc) - std::log(pc_tp1)) * (d2sat2 - d2sat1) / std::log(pc_tp2 / pc_tp1);
727 : }
728 208638 : return d2sat;
729 : }
|
