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 "PorousFlowHystereticRelativePermeabilityLiquid.h"
11 : #include "PorousFlowVanGenuchten.h"
12 :
13 : registerMooseObject("PorousFlowApp", PorousFlowHystereticRelativePermeabilityLiquid);
14 :
15 : InputParameters
16 258 : PorousFlowHystereticRelativePermeabilityLiquid::validParams()
17 : {
18 258 : InputParameters params = PorousFlowHystereticRelativePermeabilityBase::validParams();
19 774 : params.addRangeCheckedParam<Real>(
20 : "liquid_modification_range",
21 516 : 0.9,
22 : "liquid_modification_range > 0 & liquid_modification_range <= 1",
23 : "The wetting liquid relative permeability will be a cubic between S_l = "
24 : "liquid_modification_range * (1 - S_gr_Del) and S_l = 1.0 - 0.5 * S_gr_Del");
25 258 : params.addClassDescription(
26 : "PorousFlow material that computes relative permeability of the liquid phase in 1-phase or "
27 : "2-phase models that include hysteresis. You should ensure that the 'phase' for this "
28 : "Material does indeed represent the liquid phase");
29 258 : return params;
30 0 : }
31 :
32 200 : PorousFlowHystereticRelativePermeabilityLiquid::PorousFlowHystereticRelativePermeabilityLiquid(
33 200 : const InputParameters & parameters)
34 : : PorousFlowHystereticRelativePermeabilityBase(parameters),
35 198 : _liquid_phase(_phase_num),
36 198 : _liquid_modification_range(getParam<Real>("liquid_modification_range")),
37 198 : _kl_begin(_nodal_material ? declareProperty<Real>("PorousFlow_hysteresis_kr_l_begin_nodal")
38 378 : : declareProperty<Real>("PorousFlow_hysteresis_kr_l_begin_qp")),
39 198 : _klp_begin(_nodal_material ? declareProperty<Real>("PorousFlow_hysteresis_krp_l_begin_nodal")
40 378 : : declareProperty<Real>("PorousFlow_hysteresis_krp_l_begin_qp")),
41 198 : _kl_end(_nodal_material ? declareProperty<Real>("PorousFlow_hysteresis_kr_l_end_nodal")
42 378 : : declareProperty<Real>("PorousFlow_hysteresis_kr_l_end_qp")),
43 198 : _klp_end(_nodal_material ? declareProperty<Real>("PorousFlow_hysteresis_krp_l_end_nodal")
44 578 : : declareProperty<Real>("PorousFlow_hysteresis_krp_l_end_qp"))
45 : {
46 198 : }
47 :
48 : void
49 27406 : PorousFlowHystereticRelativePermeabilityLiquid::computeRelPermQp()
50 : {
51 27406 : const Real sl = _saturation[_qp][_liquid_phase];
52 :
53 27406 : if (_hys_order[_qp] == 0)
54 : {
55 14126 : _relative_permeability[_qp] = PorousFlowVanGenuchten::relativePermeabilityHys(
56 14126 : sl, _s_lr, 0.0, _s_gr_max, 1.0, _m, _liquid_modification_range, 0.0, 0.0, 0.0, 0.0);
57 14126 : _drelative_permeability_ds[_qp] = PorousFlowVanGenuchten::drelativePermeabilityHys(
58 14126 : sl, _s_lr, 0.0, _s_gr_max, 1.0, _m, _liquid_modification_range, 0.0, 0.0, 0.0, 0.0);
59 : }
60 : else
61 : {
62 : // following ternary deals with the case where the turning-point saturation occurs in the
63 : // low-saturation region (tp_sat < _s_lr). There is "no hysteresis along the extension"
64 : // according to Doughty2008, so assume that the wetting curve is the same as would occur if
65 : // the turning-point saturation occured at _s_lr
66 : const Real effective_liquid_tp =
67 13280 : (_hys_sat_tps[_qp].at(0) < _s_lr) ? _s_lr : _hys_sat_tps[_qp].at(0);
68 13280 : const Real s_gas_max = (_hys_sat_tps[_qp].at(0) < _s_lr) ? _s_gr_max : _s_gr_tp0[_qp];
69 13280 : _relative_permeability[_qp] =
70 13280 : PorousFlowVanGenuchten::relativePermeabilityHys(sl,
71 : _s_lr,
72 : s_gas_max,
73 13280 : _s_gr_max,
74 : effective_liquid_tp,
75 13280 : _m,
76 13280 : _liquid_modification_range,
77 13280 : _kl_begin[_qp],
78 13280 : _klp_begin[_qp],
79 13280 : _kl_end[_qp],
80 13280 : _klp_end[_qp]);
81 13280 : _drelative_permeability_ds[_qp] =
82 13280 : PorousFlowVanGenuchten::drelativePermeabilityHys(sl,
83 13280 : _s_lr,
84 : s_gas_max,
85 13280 : _s_gr_max,
86 : effective_liquid_tp,
87 13280 : _m,
88 13280 : _liquid_modification_range,
89 13280 : _kl_begin[_qp],
90 13280 : _klp_begin[_qp],
91 13280 : _kl_end[_qp],
92 13280 : _klp_end[_qp]);
93 : }
94 27406 : }
95 :
96 : void
97 1340 : PorousFlowHystereticRelativePermeabilityLiquid::computeTurningPoint0Info(Real tp_sat)
98 : {
99 1340 : PorousFlowHystereticRelativePermeabilityBase::computeTurningPoint0Info(tp_sat);
100 :
101 : // following ternaries deal with the case where the turning-point saturation occurs in the
102 : // low-saturation region (tp_sat < _s_lr). There is "no hysteresis along the extension"
103 : // according to Doughty2008, so assume that the wetting curve is the same as would occur if
104 : // the turning-point saturation occured at _s_lr
105 1340 : const Real effective_liquid_tp = (tp_sat < _s_lr) ? _s_lr : tp_sat;
106 1340 : const Real s_gas_max = (tp_sat < _s_lr) ? _s_gr_max : _s_gr_tp0[_qp];
107 :
108 : // compute and record the wetting-curve value and derivative at s_begin
109 1340 : const Real s_begin = _liquid_modification_range * (1.0 - s_gas_max);
110 1340 : _kl_begin[_qp] = PorousFlowVanGenuchten::relativePermeabilityHys(s_begin,
111 : _s_lr,
112 : s_gas_max,
113 1340 : _s_gr_max,
114 : effective_liquid_tp,
115 1340 : _m,
116 : _liquid_modification_range,
117 : 0.0,
118 : 0.0,
119 : 0.0,
120 : 0.0);
121 1340 : _klp_begin[_qp] = PorousFlowVanGenuchten::drelativePermeabilityHys(s_begin,
122 1340 : _s_lr,
123 : s_gas_max,
124 1340 : _s_gr_max,
125 : effective_liquid_tp,
126 1340 : _m,
127 1340 : _liquid_modification_range,
128 : 0.0,
129 : 0.0,
130 : 0.0,
131 : 0.0);
132 : // compute and record the drying curve information at s_end (the drying curve is used because
133 : // s_end = 1.0 - 0.5 * s_gas_max in PorousFlowVanGenuchten::relativePermeabilityHys)
134 1340 : const Real s_end = 1.0 - 0.5 * s_gas_max;
135 1340 : _kl_end[_qp] = PorousFlowVanGenuchten::relativePermeabilityHys(s_end,
136 1340 : _s_lr,
137 : s_gas_max,
138 1340 : _s_gr_max,
139 : effective_liquid_tp,
140 1340 : _m,
141 1340 : _liquid_modification_range,
142 : 0.0,
143 : 0.0,
144 : 0.0,
145 : 0.0);
146 2680 : _klp_end[_qp] = PorousFlowVanGenuchten::drelativePermeabilityHys(s_end,
147 1340 : _s_lr,
148 : s_gas_max,
149 1340 : _s_gr_max,
150 : effective_liquid_tp,
151 1340 : _m,
152 1340 : _liquid_modification_range,
153 : 0.0,
154 : 0.0,
155 : 0.0,
156 : 0.0);
157 1340 : }
|
