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 4619 : PorousFlowHystereticCapillaryPressure::validParams()
17 : {
18 4619 : InputParameters params = PorousFlowVariableBase::validParams();
19 9238 : 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 9238 : 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 9238 : 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 9238 : 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 9238 : 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 13857 : params.addRangeCheckedParam<Real>(
43 : "S_lr",
44 9238 : 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 9238 : 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 13857 : params.addRangeCheckedParam<Real>(
57 : "Pc_max",
58 9238 : 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 9238 : MooseEnum low_ext_enum("none quadratic exponential", "exponential");
64 9238 : 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 13857 : params.addRangeCheckedParam<Real>(
77 : "high_ratio",
78 9238 : 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 9238 : MooseEnum high_ext_enum("none power", "power");
83 9238 : 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 4619 : 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 4619 : return params;
97 4619 : }
98 :
99 3558 : PorousFlowHystereticCapillaryPressure::PorousFlowHystereticCapillaryPressure(
100 3558 : const InputParameters & parameters)
101 : : PorousFlowVariableBase(parameters),
102 3558 : _alpha_d(getParam<Real>("alpha_d")),
103 7116 : _alpha_w(getParam<Real>("alpha_w")),
104 7116 : _n_d(getParam<Real>("n_d")),
105 7116 : _n_w(getParam<Real>("n_w")),
106 7116 : _s_l_min(getParam<Real>("S_l_min")),
107 7116 : _s_lr(getParam<Real>("S_lr")),
108 7116 : _s_gr_max(getParam<Real>("S_gr_max")),
109 7116 : _pc_max(getParam<Real>("Pc_max")),
110 7116 : _high_ratio(getParam<Real>("high_ratio")),
111 3558 : _low_ext_type(
112 7116 : getParam<MooseEnum>("low_extension_type")
113 : .getEnum<PorousFlowVanGenuchten::LowCapillaryPressureExtension::ExtensionStrategy>()),
114 3558 : _s_low_d(PorousFlowVanGenuchten::saturationHys(_pc_max, _s_l_min, 0.0, _alpha_d, _n_d)),
115 3558 : _dpc_low_d(
116 3558 : PorousFlowVanGenuchten::dcapillaryPressureHys(_s_low_d, _s_l_min, 0.0, _alpha_d, _n_d)),
117 3558 : _low_ext_d(_low_ext_type, _s_low_d, _pc_max, _dpc_low_d),
118 3558 : _s_low_w(PorousFlowVanGenuchten::saturationHys(_pc_max, _s_l_min, _s_gr_max, _alpha_w, _n_w)),
119 7116 : _dpc_low_w(PorousFlowVanGenuchten::dcapillaryPressureHys(
120 3558 : _s_low_w, _s_l_min, _s_gr_max, _alpha_w, _n_w)),
121 3558 : _low_ext_w(_low_ext_type, _s_low_w, _pc_max, _dpc_low_w),
122 3558 : _high_ext_type(
123 3558 : getParam<MooseEnum>("high_extension_type")
124 : .getEnum<PorousFlowVanGenuchten::HighCapillaryPressureExtension::ExtensionStrategy>()),
125 3558 : _s_high(_high_ratio * (1 - _s_gr_max)),
126 3558 : _pc_high(
127 3558 : PorousFlowVanGenuchten::capillaryPressureHys(_s_high, _s_l_min, _s_gr_max, _alpha_w, _n_w)),
128 7116 : _dpc_high(PorousFlowVanGenuchten::dcapillaryPressureHys(
129 3558 : _s_high, _s_l_min, _s_gr_max, _alpha_w, _n_w)),
130 3558 : _high_ext(_high_ext_type, _s_high, _pc_high, _dpc_high),
131 3558 : _hys_order(_nodal_material ? getMaterialProperty<unsigned>("PorousFlow_hysteresis_order_nodal")
132 10210 : : getMaterialProperty<unsigned>("PorousFlow_hysteresis_order_qp")),
133 7116 : _hys_order_old(_nodal_material
134 3558 : ? getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_nodal")
135 10210 : : getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_qp")),
136 3558 : _hys_sat_tps(
137 3558 : _nodal_material
138 3558 : ? getMaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
139 : "PorousFlow_hysteresis_saturation_tps_nodal")
140 10210 : : getMaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
141 : "PorousFlow_hysteresis_saturation_tps_qp")),
142 7116 : _pc_older(_nodal_material
143 3558 : ? getMaterialPropertyOlder<Real>("PorousFlow_hysteretic_capillary_pressure_nodal")
144 10210 : : getMaterialPropertyOlder<Real>("PorousFlow_hysteretic_capillary_pressure_qp")),
145 7116 : _pc_tps(_nodal_material
146 3558 : ? declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
147 : "PorousFlow_hysteresis_Pc_tps_nodal")
148 6884 : : declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
149 : "PorousFlow_hysteresis_Pc_tps_qp")),
150 7116 : _s_d_tps(_nodal_material
151 3558 : ? declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
152 : "PorousFlow_hysteresis_s_d_tps_nodal")
153 6884 : : declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
154 : "PorousFlow_hysteresis_s_d_tps_qp")),
155 7116 : _s_gr_tps(_nodal_material
156 3558 : ? declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
157 : "PorousFlow_hysteresis_s_gr_tps_nodal")
158 6884 : : declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
159 : "PorousFlow_hysteresis_s_gr_tps_qp")),
160 3558 : _w_low_ext_tps(
161 3558 : _nodal_material
162 3558 : ? declareProperty<std::array<PorousFlowVanGenuchten::LowCapillaryPressureExtension,
163 232 : PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
164 : "PorousFlow_hysteresis_w_low_ext_tps_nodal")
165 : : declareProperty<std::array<PorousFlowVanGenuchten::LowCapillaryPressureExtension,
166 6884 : PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
167 : "PorousFlow_hysteresis_w_low_ext_tps_qp")),
168 3558 : _w_high_ext_tps(
169 3558 : _nodal_material
170 3558 : ? declareProperty<std::array<PorousFlowVanGenuchten::HighCapillaryPressureExtension,
171 232 : PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
172 : "PorousFlow_hysteresis_w_high_ext_tps_nodal")
173 : : declareProperty<std::array<PorousFlowVanGenuchten::HighCapillaryPressureExtension,
174 6884 : PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
175 : "PorousFlow_hysteresis_w_high_ext_tps_qp")),
176 7116 : _s_w_tps(_nodal_material
177 3558 : ? declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
178 : "PorousFlow_hysteresis_s_w_tps_nodal")
179 6884 : : declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
180 3558 : "PorousFlow_hysteresis_s_w_tps_qp"))
181 : {
182 3558 : if (_s_gr_max >= 1 - _s_l_min)
183 2 : paramError("S_gr_max", "Must be less than 1 - S_l_min");
184 3556 : 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 3554 : if (_s_lr <= _s_l_min)
189 2 : mooseWarning("S_lr should usually be greater than S_l_min");
190 3552 : }
191 :
192 : void
193 213024 : PorousFlowHystereticCapillaryPressure::initQpStatefulProperties()
194 : {
195 213024 : PorousFlowVariableBase::initQpStatefulProperties();
196 : Real tp_sat;
197 : Real tp_pc;
198 213024 : if (_hys_order[_qp] >= 1)
199 : {
200 157836 : tp_sat = _hys_sat_tps[_qp].at(0);
201 315672 : tp_pc = PorousFlowVanGenuchten::capillaryPressureHys(
202 157836 : tp_sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
203 157836 : computeTurningPointInfo(0, tp_sat, tp_pc);
204 : }
205 213024 : if (_hys_order[_qp] >= 2)
206 : {
207 79820 : tp_sat = _hys_sat_tps[_qp].at(1);
208 79820 : tp_pc = firstOrderWettingPc(tp_sat);
209 79820 : computeTurningPointInfo(1, tp_sat, tp_pc);
210 : }
211 213024 : if (_hys_order[_qp] >= 3)
212 : {
213 29900 : tp_sat = _hys_sat_tps[_qp].at(2);
214 29900 : tp_pc = secondOrderDryingPc(tp_sat);
215 29900 : computeTurningPointInfo(2, tp_sat, tp_pc);
216 : }
217 213024 : }
218 :
219 : void
220 573258 : PorousFlowHystereticCapillaryPressure::computeQpProperties()
221 : {
222 : // size stuff correctly and prepare the derivative matrices with zeroes
223 573258 : PorousFlowVariableBase::computeQpProperties();
224 :
225 573258 : if (_hys_order[_qp] != _hys_order_old[_qp] && _hys_order[_qp] > 0)
226 : {
227 9104 : const unsigned tp_num = _hys_order[_qp] - 1;
228 9104 : computeTurningPointInfo(tp_num, _hys_sat_tps[_qp].at(tp_num), _pc_older[_qp]);
229 : }
230 573258 : }
231 :
232 : Real
233 276660 : PorousFlowHystereticCapillaryPressure::landSat(Real slDel) const
234 : {
235 276660 : const Real a = 1.0 / _s_gr_max - 1.0 / (1.0 - _s_lr);
236 276660 : return (1.0 - slDel) / (1.0 + a * (1.0 - slDel));
237 : }
238 :
239 : void
240 276660 : PorousFlowHystereticCapillaryPressure::computeTurningPointInfo(unsigned tp_num,
241 : Real tp_sat,
242 : Real tp_pc)
243 : {
244 276660 : _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 276660 : const Real pc_d_tps = PorousFlowVanGenuchten::capillaryPressureHys(
249 276660 : tp_sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
250 : // saturation on the drying curve, at tp_pc
251 276660 : _s_d_tps[_qp].at(tp_num) =
252 276660 : 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 276660 : _s_gr_tps[_qp].at(tp_num) = landSat(tp_sat);
257 : // the low extension of the wetting curve defined by _s_gr_tps
258 276660 : const Real s_w_low_ext = PorousFlowVanGenuchten::saturationHys(
259 276660 : _pc_max, _s_l_min, _s_gr_tps[_qp].at(tp_num), _alpha_w, _n_w);
260 276660 : const Real dpc_w_low_ext = PorousFlowVanGenuchten::dcapillaryPressureHys(
261 276660 : s_w_low_ext, _s_l_min, _s_gr_tps[_qp].at(tp_num), _alpha_w, _n_w);
262 276660 : _w_low_ext_tps[_qp].at(tp_num) = PorousFlowVanGenuchten::LowCapillaryPressureExtension(
263 276660 : _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 276660 : (_high_ext_type == PorousFlowVanGenuchten::HighCapillaryPressureExtension::NONE)
267 276660 : ? 1.0 - _s_gr_tps[_qp].at(tp_num)
268 166160 : : _high_ratio *
269 166160 : (1.0 -
270 166160 : _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 276660 : PorousFlowVanGenuchten::capillaryPressureHys(s_w_high_ext,
274 276660 : _s_l_min,
275 276660 : _s_gr_tps[_qp].at(tp_num),
276 276660 : _alpha_w,
277 276660 : _n_w,
278 : _w_low_ext_tps[_qp].at(tp_num));
279 : const Real dpc_w_high_ext =
280 276660 : PorousFlowVanGenuchten::dcapillaryPressureHys(s_w_high_ext,
281 276660 : _s_l_min,
282 276660 : _s_gr_tps[_qp].at(tp_num),
283 276660 : _alpha_w,
284 276660 : _n_w,
285 276660 : _w_low_ext_tps[_qp].at(tp_num));
286 276660 : _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 553320 : _s_w_tps[_qp].at(tp_num) = PorousFlowVanGenuchten::saturationHys(pc_d_tps,
291 276660 : _s_l_min,
292 276660 : _s_gr_tps[_qp].at(tp_num),
293 276660 : _alpha_w,
294 276660 : _n_w,
295 276660 : _w_low_ext_tps[_qp].at(tp_num),
296 : _w_high_ext_tps[_qp].at(tp_num));
297 276660 : }
298 :
299 : Real
300 1322184 : PorousFlowHystereticCapillaryPressure::capillaryPressureQp(Real sat) const
301 : {
302 : Real pc = 0.0;
303 1322184 : if (_hys_order[_qp] == 0) // on primary drying curve
304 313136 : pc = PorousFlowVanGenuchten::capillaryPressureHys(
305 313136 : sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
306 1009048 : else if (_hys_order[_qp] == 1) // first-order wetting
307 501584 : pc = firstOrderWettingPc(sat);
308 507464 : else if (_hys_order[_qp] == 2) // second-order drying
309 317904 : pc = secondOrderDryingPc(sat);
310 : else // third order drying and wetting
311 : {
312 189560 : const Real tp1 = _hys_sat_tps[_qp].at(1);
313 189560 : const Real tp2 = _hys_sat_tps[_qp].at(2);
314 189560 : const Real pc1 = firstOrderWettingPc(sat);
315 189560 : 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 189560 : if (sat < tp2)
318 : pc = pc2;
319 189560 : else if (sat > tp1)
320 : pc = pc1;
321 189464 : else if (pc1 >= pc2)
322 : pc = pc1;
323 189464 : else if (pc1 <= 0.0 || pc2 <= 0.0)
324 : pc = 0.0;
325 : else
326 157730 : pc = std::exp(std::log(pc1) + (sat - tp1) * (std::log(pc2) - std::log(pc1)) / (tp2 - tp1));
327 : }
328 1322184 : return pc;
329 : }
330 :
331 : Real
332 262872 : PorousFlowHystereticCapillaryPressure::dcapillaryPressureQp(Real sat) const
333 : {
334 : Real dpc = 0.0;
335 262872 : if (_hys_order[_qp] == 0) // on primary drying curve
336 85368 : dpc = PorousFlowVanGenuchten::dcapillaryPressureHys(
337 85368 : sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
338 177504 : else if (_hys_order[_qp] == 1) // first-order wetting
339 79980 : dpc = dfirstOrderWettingPc(sat);
340 97524 : else if (_hys_order[_qp] == 2) // second-order drying
341 60624 : dpc = dsecondOrderDryingPc(sat);
342 : else // third order drying and wetting
343 : {
344 36900 : const Real tp1 = _hys_sat_tps[_qp].at(1);
345 36900 : const Real tp2 = _hys_sat_tps[_qp].at(2);
346 36900 : const Real pc1 = firstOrderWettingPc(sat);
347 36900 : 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 36900 : if (sat < tp2)
350 0 : dpc = dsecondOrderDryingPc(sat);
351 36900 : else if (sat > tp1)
352 0 : dpc = dfirstOrderWettingPc(sat);
353 36900 : else if (pc1 >= pc2)
354 0 : dpc = dfirstOrderWettingPc(sat);
355 36900 : else if (pc1 <= 0.0 || pc2 <= 0.0)
356 : dpc = 0.0;
357 : else
358 : {
359 : const Real pc =
360 31611 : std::exp(std::log(pc1) + (sat - tp1) * (std::log(pc2) - std::log(pc1)) / (tp2 - tp1));
361 31611 : const Real dpc1 = dfirstOrderWettingPc(sat);
362 31611 : const Real dpc2 = dsecondOrderDryingPc(sat);
363 31611 : dpc = pc * (dpc1 / pc1 + (std::log(pc2) - std::log(pc1)) / (tp2 - tp1) +
364 31611 : (sat - tp1) * (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1));
365 : }
366 : }
367 262872 : return dpc;
368 : }
369 :
370 : Real
371 90846 : PorousFlowHystereticCapillaryPressure::d2capillaryPressureQp(Real sat) const
372 : {
373 : Real d2pc = 0.0;
374 90846 : if (_hys_order[_qp] == 0) // on primary drying curve
375 30384 : d2pc = PorousFlowVanGenuchten::d2capillaryPressureHys(
376 30384 : sat, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
377 60462 : else if (_hys_order[_qp] == 1) // first-order wetting
378 27690 : d2pc = d2firstOrderWettingPc(sat);
379 32772 : else if (_hys_order[_qp] == 2) // second-order drying
380 20472 : d2pc = d2secondOrderDryingPc(sat);
381 : else // third order drying and wetting
382 : {
383 12300 : const Real tp1 = _hys_sat_tps[_qp].at(1);
384 12300 : const Real tp2 = _hys_sat_tps[_qp].at(2);
385 12300 : const Real pc1 = firstOrderWettingPc(sat);
386 12300 : 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 12300 : if (sat < tp2)
389 0 : d2pc = d2secondOrderDryingPc(sat);
390 12300 : else if (sat > tp1)
391 0 : d2pc = d2firstOrderWettingPc(sat);
392 12300 : else if (pc1 >= pc2)
393 0 : d2pc = d2firstOrderWettingPc(sat);
394 12300 : else if (pc1 <= 0.0 || pc2 <= 0.0)
395 : d2pc = 0.0;
396 : else
397 : {
398 : const Real pc =
399 10537 : std::exp(std::log(pc1) + (sat - tp1) * (std::log(pc2) - std::log(pc1)) / (tp2 - tp1));
400 10537 : const Real dpc1 = dfirstOrderWettingPc(sat);
401 10537 : const Real dpc2 = dsecondOrderDryingPc(sat);
402 10537 : const Real dpc = pc * (dpc1 / pc1 + (std::log(pc2) - std::log(pc1)) / (tp2 - tp1) +
403 10537 : (sat - tp1) * (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1));
404 10537 : const Real d2pc1 = d2firstOrderWettingPc(sat);
405 10537 : const Real d2pc2 = d2secondOrderDryingPc(sat);
406 10537 : d2pc =
407 10537 : dpc * (dpc1 / pc1 + (std::log(pc2) - std::log(pc1)) / (tp2 - tp1) +
408 : (sat - tp1) * (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1)) +
409 10537 : pc *
410 10537 : (d2pc1 / pc1 - dpc1 * dpc1 / pc1 / pc1 + (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1) +
411 10537 : (dpc2 / pc2 - dpc1 / pc1) / (tp2 - tp1) +
412 10537 : (sat - tp1) *
413 10537 : (d2pc2 / pc2 - dpc2 * dpc2 / pc2 / pc2 - d2pc1 / pc1 + dpc1 * dpc1 / pc1 / pc1) /
414 : (tp2 - tp1));
415 : }
416 : }
417 90846 : return d2pc;
418 : }
419 :
420 : Real
421 820164 : 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 820164 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
427 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
428 820164 : ? 1.0
429 335900 : : 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 820164 : const Real sat_to_use = _s_w_tps[_qp].at(0) + (max_s - _s_w_tps[_qp].at(0)) *
434 820164 : (sat - _hys_sat_tps[_qp].at(0)) /
435 820164 : (max_s - _hys_sat_tps[_qp].at(0));
436 820164 : return PorousFlowVanGenuchten::capillaryPressureHys(sat_to_use,
437 820164 : _s_l_min,
438 820164 : _s_gr_tps[_qp].at(0),
439 820164 : _alpha_w,
440 820164 : _n_w,
441 820164 : _w_low_ext_tps[_qp].at(0),
442 820164 : _w_high_ext_tps[_qp].at(0));
443 : }
444 :
445 : Real
446 122128 : PorousFlowHystereticCapillaryPressure::dfirstOrderWettingPc(Real sat) const
447 : {
448 122128 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
449 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
450 122128 : ? 1.0
451 33948 : : 1.0 - _s_gr_tps[_qp].at(0);
452 122128 : const Real sat_to_use = _s_w_tps[_qp].at(0) + (max_s - _s_w_tps[_qp].at(0)) *
453 122128 : (sat - _hys_sat_tps[_qp].at(0)) /
454 122128 : (max_s - _hys_sat_tps[_qp].at(0));
455 122128 : const Real dsat_to_use = (max_s - _s_w_tps[_qp].at(0)) / (max_s - _hys_sat_tps[_qp].at(0));
456 122128 : return PorousFlowVanGenuchten::dcapillaryPressureHys(sat_to_use,
457 122128 : _s_l_min,
458 122128 : _s_gr_tps[_qp].at(0),
459 122128 : _alpha_w,
460 122128 : _n_w,
461 122128 : _w_low_ext_tps[_qp].at(0),
462 : _w_high_ext_tps[_qp].at(0)) *
463 122128 : dsat_to_use;
464 : }
465 :
466 : Real
467 38227 : PorousFlowHystereticCapillaryPressure::d2firstOrderWettingPc(Real sat) const
468 : {
469 38227 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
470 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
471 38227 : ? 1.0
472 10537 : : 1.0 - _s_gr_tps[_qp].at(0);
473 38227 : const Real sat_to_use = _s_w_tps[_qp].at(0) + (max_s - _s_w_tps[_qp].at(0)) *
474 38227 : (sat - _hys_sat_tps[_qp].at(0)) /
475 38227 : (max_s - _hys_sat_tps[_qp].at(0));
476 38227 : const Real dsat_to_use = (max_s - _s_w_tps[_qp].at(0)) / (max_s - _hys_sat_tps[_qp].at(0));
477 38227 : return PorousFlowVanGenuchten::d2capillaryPressureHys(sat_to_use,
478 38227 : _s_l_min,
479 38227 : _s_gr_tps[_qp].at(0),
480 38227 : _alpha_w,
481 38227 : _n_w,
482 38227 : _w_low_ext_tps[_qp].at(0),
483 38227 : _w_high_ext_tps[_qp].at(0)) *
484 38227 : dsat_to_use * dsat_to_use;
485 : }
486 :
487 : Real
488 586564 : 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 586564 : const Real tp0 = _hys_sat_tps[_qp].at(0);
494 586564 : const Real tp1 = _hys_sat_tps[_qp].at(1);
495 586564 : const Real s1 = _s_d_tps[_qp].at(1);
496 : const Real sat_to_use =
497 586564 : (sat >= tp0) ? tp0 + (sat - tp0) * (s1 - tp0) / (tp1 - tp0)
498 : : sat; // final case can occur just at transition from 2nd to 0th order
499 1173128 : return PorousFlowVanGenuchten::capillaryPressureHys(
500 586564 : sat_to_use, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
501 : }
502 :
503 : Real
504 102772 : PorousFlowHystereticCapillaryPressure::dsecondOrderDryingPc(Real sat) const
505 : {
506 102772 : const Real tp0 = _hys_sat_tps[_qp].at(0);
507 102772 : const Real tp1 = _hys_sat_tps[_qp].at(1);
508 102772 : const Real s1 = _s_d_tps[_qp].at(1);
509 102772 : const Real sat_to_use = (sat >= tp0) ? tp0 + (sat - tp0) * (s1 - tp0) / (tp1 - tp0) : sat;
510 102772 : const Real dsat_to_use = (sat >= tp0) ? (s1 - tp0) / (tp1 - tp0) : 1.0;
511 205544 : return PorousFlowVanGenuchten::dcapillaryPressureHys(
512 102772 : sat_to_use, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d) *
513 102772 : dsat_to_use;
514 : }
515 :
516 : Real
517 31009 : PorousFlowHystereticCapillaryPressure::d2secondOrderDryingPc(Real sat) const
518 : {
519 31009 : const Real tp0 = _hys_sat_tps[_qp].at(0);
520 31009 : const Real tp1 = _hys_sat_tps[_qp].at(1);
521 31009 : const Real s1 = _s_d_tps[_qp].at(1);
522 31009 : const Real sat_to_use = (sat >= tp0) ? tp0 + (sat - tp0) * (s1 - tp0) / (tp1 - tp0) : sat;
523 31009 : const Real dsat_to_use = (sat >= tp0) ? (s1 - tp0) / (tp1 - tp0) : 1.0;
524 62018 : return PorousFlowVanGenuchten::d2capillaryPressureHys(
525 31009 : sat_to_use, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d) *
526 31009 : dsat_to_use * dsat_to_use;
527 : }
528 :
529 : Real
530 241020 : PorousFlowHystereticCapillaryPressure::firstOrderWettingSat(Real pc) const
531 : {
532 : // this is inverse of firstOrderWettingPc: see that method for comments
533 241020 : const Real sat_to_use = PorousFlowVanGenuchten::saturationHys(pc,
534 241020 : _s_l_min,
535 241020 : _s_gr_tps[_qp].at(0),
536 241020 : _alpha_w,
537 241020 : _n_w,
538 241020 : _w_low_ext_tps[_qp].at(0),
539 241020 : _w_high_ext_tps[_qp].at(0));
540 241020 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
541 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
542 241020 : ? 1.0
543 73800 : : 1.0 - _s_gr_tps[_qp].at(0);
544 241020 : if (sat_to_use > max_s) // this occurs when using no high extension and pc = 0
545 : return max_s;
546 220069 : return (sat_to_use - _s_w_tps[_qp].at(0)) * (max_s - _hys_sat_tps[_qp].at(0)) /
547 220069 : (max_s - _s_w_tps[_qp].at(0)) +
548 220069 : _hys_sat_tps[_qp].at(0);
549 : }
550 :
551 : Real
552 166656 : PorousFlowHystereticCapillaryPressure::dfirstOrderWettingSat(Real pc) const
553 : {
554 166656 : const Real dsat_to_use = PorousFlowVanGenuchten::dsaturationHys(pc,
555 166656 : _s_l_min,
556 166656 : _s_gr_tps[_qp].at(0),
557 166656 : _alpha_w,
558 166656 : _n_w,
559 166656 : _w_low_ext_tps[_qp].at(0),
560 166656 : _w_high_ext_tps[_qp].at(0));
561 166656 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
562 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
563 166656 : ? 1.0
564 41000 : : 1.0 - _s_gr_tps[_qp].at(0);
565 166656 : if (pc <= 0.0 && _w_high_ext_tps[_qp].at(0).strategy ==
566 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::NONE)
567 : return 0.0;
568 156119 : 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 56664 : PorousFlowHystereticCapillaryPressure::d2firstOrderWettingSat(Real pc) const
573 : {
574 56664 : const Real d2sat_to_use = PorousFlowVanGenuchten::d2saturationHys(pc,
575 56664 : _s_l_min,
576 56664 : _s_gr_tps[_qp].at(0),
577 56664 : _alpha_w,
578 56664 : _n_w,
579 56664 : _w_low_ext_tps[_qp].at(0),
580 56664 : _w_high_ext_tps[_qp].at(0));
581 56664 : const Real max_s = (_w_high_ext_tps[_qp].at(0).strategy ==
582 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::POWER)
583 56664 : ? 1.0
584 12300 : : 1.0 - _s_gr_tps[_qp].at(0);
585 56664 : if (pc <= 0.0 && _w_high_ext_tps[_qp].at(0).strategy ==
586 : PorousFlowVanGenuchten::HighCapillaryPressureExtension::NONE)
587 : return 0.0;
588 52277 : 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 191572 : 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 191572 : PorousFlowVanGenuchten::saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
597 191572 : const Real tp0 = _hys_sat_tps[_qp].at(0);
598 191572 : const Real tp1 = _hys_sat_tps[_qp].at(1);
599 191572 : const Real s1 = _s_d_tps[_qp].at(1);
600 191572 : return (sat_to_use >= tp0) ? (sat_to_use - tp0) * (tp1 - tp0) / (s1 - tp0) + tp0 : sat_to_use;
601 : }
602 :
603 : Real
604 122136 : PorousFlowHystereticCapillaryPressure::dsecondOrderDryingSat(Real pc) const
605 : {
606 : const Real sat_to_use =
607 122136 : PorousFlowVanGenuchten::saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
608 : const Real dsat_to_use =
609 122136 : PorousFlowVanGenuchten::dsaturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
610 122136 : const Real tp0 = _hys_sat_tps[_qp].at(0);
611 122136 : const Real tp1 = _hys_sat_tps[_qp].at(1);
612 122136 : const Real s1 = _s_d_tps[_qp].at(1);
613 122136 : return (sat_to_use >= tp0) ? dsat_to_use * (tp1 - tp0) / (s1 - tp0) : dsat_to_use;
614 : }
615 :
616 : Real
617 37608 : PorousFlowHystereticCapillaryPressure::d2secondOrderDryingSat(Real pc) const
618 : {
619 : const Real sat_to_use =
620 37608 : PorousFlowVanGenuchten::saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
621 : const Real d2sat_to_use =
622 37608 : PorousFlowVanGenuchten::d2saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
623 37608 : const Real tp0 = _hys_sat_tps[_qp].at(0);
624 37608 : const Real tp1 = _hys_sat_tps[_qp].at(1);
625 37608 : const Real s1 = _s_d_tps[_qp].at(1);
626 37608 : return (sat_to_use >= tp0) ? d2sat_to_use * (tp1 - tp0) / (s1 - tp0) : d2sat_to_use;
627 : }
628 :
629 : Real
630 411484 : PorousFlowHystereticCapillaryPressure::liquidSaturationQp(Real pc) const
631 : {
632 : Real sat = 0.0;
633 411484 : if (_hys_order[_qp] == 0) // on primary drying curve
634 141016 : sat = PorousFlowVanGenuchten::saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
635 270468 : else if (_hys_order[_qp] == 1) // first-order wetting
636 135080 : sat = firstOrderWettingSat(pc);
637 135388 : else if (_hys_order[_qp] == 2) // second-order drying
638 85632 : 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 49756 : const Real pc_tp1 = _pc_tps[_qp].at(1); // pc on first-order wetting at TP_1
647 49756 : const Real pc_tp2 = _pc_tps[_qp].at(2); // pc on second-order drying at TP_2
648 49756 : const Real sat1 = firstOrderWettingSat(pc);
649 49756 : 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 49756 : 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 33308 : else if (pc > pc_tp2)
657 : sat = sat2;
658 33308 : else if (pc < pc_tp1)
659 : sat = sat1;
660 : else
661 32948 : sat = sat1 + (std::log(pc) - std::log(pc_tp1)) * (sat2 - sat1) / std::log(pc_tp2 / pc_tp1);
662 : }
663 411484 : return sat;
664 : }
665 :
666 : Real
667 350850 : PorousFlowHystereticCapillaryPressure::dliquidSaturationQp(Real pc) const
668 : {
669 : Real dsat = 0.0;
670 350850 : if (_hys_order[_qp] == 0) // on primary drying curve
671 132870 : dsat = PorousFlowVanGenuchten::dsaturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
672 217980 : else if (_hys_order[_qp] == 1) // first-order wetting
673 110472 : dsat = dfirstOrderWettingSat(pc);
674 107508 : else if (_hys_order[_qp] == 2) // second-order drying
675 65952 : dsat = dsecondOrderDryingSat(pc);
676 : else // third order drying and wetting
677 : {
678 41556 : const Real pc_tp1 = _pc_tps[_qp].at(1); // pc on first-order wetting at TP_1
679 41556 : const Real pc_tp2 = _pc_tps[_qp].at(2); // pc on second-order drying at TP_2
680 41556 : const Real sat1 = firstOrderWettingSat(pc);
681 41556 : const Real sat2 = secondOrderDryingSat(pc);
682 41556 : const Real dsat1 = dfirstOrderWettingSat(pc);
683 41556 : const Real dsat2 = dsecondOrderDryingSat(pc);
684 41556 : if (pc_tp1 <= 0.0 || pc <= 0.0 || pc_tp2 <= 0.0)
685 : dsat = dsat2;
686 29208 : else if (pc > pc_tp2)
687 : dsat = dsat2;
688 29208 : else if (pc < pc_tp1)
689 : dsat = dsat1;
690 : else
691 28848 : dsat = dsat1 + 1.0 / pc * (sat2 - sat1) / std::log(pc_tp2 / pc_tp1) +
692 28848 : (std::log(pc) - std::log(pc_tp1)) * (dsat2 - dsat1) / std::log(pc_tp2 / pc_tp1);
693 : }
694 350850 : return dsat;
695 : }
696 :
697 : Real
698 140946 : PorousFlowHystereticCapillaryPressure::d2liquidSaturationQp(Real pc) const
699 : {
700 : Real d2sat = 0.0;
701 140946 : if (_hys_order[_qp] == 0) // on primary drying curve
702 61302 : d2sat = PorousFlowVanGenuchten::d2saturationHys(pc, _s_l_min, 0.0, _alpha_d, _n_d, _low_ext_d);
703 79644 : else if (_hys_order[_qp] == 1) // first-order wetting
704 42036 : d2sat = d2firstOrderWettingSat(pc);
705 37608 : else if (_hys_order[_qp] == 2) // second-order drying
706 22980 : d2sat = d2secondOrderDryingSat(pc);
707 : else // third order drying and wetting
708 : {
709 14628 : const Real pc_tp1 = _pc_tps[_qp].at(1); // pc on first-order wetting at TP_1
710 14628 : const Real pc_tp2 = _pc_tps[_qp].at(2); // pc on second-order drying at TP_2
711 14628 : const Real sat1 = firstOrderWettingSat(pc);
712 14628 : const Real sat2 = secondOrderDryingSat(pc);
713 14628 : const Real dsat1 = dfirstOrderWettingSat(pc);
714 14628 : const Real dsat2 = dsecondOrderDryingSat(pc);
715 14628 : const Real d2sat1 = d2firstOrderWettingSat(pc);
716 14628 : const Real d2sat2 = d2secondOrderDryingSat(pc);
717 14628 : if (pc_tp1 <= 0.0 || pc <= 0.0 || pc_tp2 <= 0.0)
718 : d2sat = d2sat2;
719 10504 : else if (pc > pc_tp2)
720 : d2sat = d2sat2;
721 10504 : else if (pc < pc_tp1)
722 : d2sat = d2sat1;
723 : else
724 10324 : d2sat = d2sat1 - 1.0 / pc / pc * (sat2 - sat1) / std::log(pc_tp2 / pc_tp1) +
725 10324 : 2.0 / pc * (dsat2 - dsat1) / std::log(pc_tp2 / pc_tp1) +
726 10324 : (std::log(pc) - std::log(pc_tp1)) * (d2sat2 - d2sat1) / std::log(pc_tp2 / pc_tp1);
727 : }
728 140946 : return d2sat;
729 : }
|
