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 "PorousFlowUnsaturated.h"
11 :
12 : #include "FEProblem.h"
13 : #include "Conversion.h"
14 : #include "libmesh/string_to_enum.h"
15 :
16 : registerMooseAction("PorousFlowApp", PorousFlowUnsaturated, "add_user_object");
17 :
18 : registerMooseAction("PorousFlowApp", PorousFlowUnsaturated, "add_kernel");
19 :
20 : registerMooseAction("PorousFlowApp", PorousFlowUnsaturated, "add_material");
21 :
22 : registerMooseAction("PorousFlowApp", PorousFlowUnsaturated, "add_aux_variable");
23 :
24 : registerMooseAction("PorousFlowApp", PorousFlowUnsaturated, "add_aux_kernel");
25 :
26 : InputParameters
27 784 : PorousFlowUnsaturated::validParams()
28 : {
29 784 : InputParameters params = PorousFlowSinglePhaseBase::validParams();
30 1568 : params.addParam<bool>("add_saturation_aux", true, "Add an AuxVariable that records saturation");
31 2352 : params.addRangeCheckedParam<Real>(
32 : "van_genuchten_alpha",
33 1568 : 1.0E-6,
34 : "van_genuchten_alpha > 0.0",
35 : "Van Genuchten alpha parameter used to determine saturation from porepressure");
36 2352 : params.addRangeCheckedParam<Real>(
37 : "van_genuchten_m",
38 1568 : 0.6,
39 : "van_genuchten_m > 0 & van_genuchten_m < 1",
40 : "Van Genuchten m parameter used to determine saturation from porepressure");
41 1568 : MooseEnum relperm_type_choice("FLAC Corey", "FLAC");
42 1568 : params.addParam<MooseEnum>("relative_permeability_type",
43 : relperm_type_choice,
44 : "Type of relative-permeability function. FLAC relperm = (1+m)S^m - "
45 : "mS^(1+m). Corey relperm = S^m. m is the exponent. Here S = "
46 : "(saturation - residual)/(1 - residual)");
47 2352 : params.addRangeCheckedParam<Real>("relative_permeability_exponent",
48 1568 : 3.0,
49 : "relative_permeability_exponent>=0",
50 : "Relative permeability exponent");
51 2352 : params.addRangeCheckedParam<Real>(
52 : "residual_saturation",
53 1568 : 0.0,
54 : "residual_saturation>=0.0 & residual_saturation<1.0",
55 : "Residual saturation to use in the relative permeability expression");
56 784 : params.addClassDescription("Adds Kernels and fluid-property Materials necessary to simulate a "
57 : "single-phase saturated-unsaturated flow problem. The saturation is "
58 : "computed using van Genuchten's expression. No Kernels for diffusion "
59 : "and dispersion of fluid components are added. To run a simulation "
60 : "you will also need to provide various other Materials for each mesh "
61 : "block, depending on your simulation type, viz: permeability, "
62 : "porosity, elasticity tensor, strain calculator, stress calculator, "
63 : "matrix internal energy, thermal conductivity, diffusivity");
64 784 : return params;
65 784 : }
66 :
67 784 : PorousFlowUnsaturated::PorousFlowUnsaturated(const InputParameters & params)
68 : : PorousFlowSinglePhaseBase(params),
69 782 : _add_saturation_aux(getParam<bool>("add_saturation_aux")),
70 1564 : _van_genuchten_alpha(getParam<Real>("van_genuchten_alpha")),
71 1564 : _van_genuchten_m(getParam<Real>("van_genuchten_m")),
72 782 : _relperm_type(
73 782 : getParam<MooseEnum>("relative_permeability_type").getEnum<RelpermTypeChoiceEnum>()),
74 1564 : _relative_permeability_exponent(getParam<Real>("relative_permeability_exponent")),
75 1564 : _s_res(getParam<Real>("residual_saturation")),
76 1566 : _capillary_pressure_name("PorousFlowUnsaturated_CapillaryPressureVG")
77 : {
78 782 : if (_stabilization == StabilizationEnum::None)
79 2 : paramError("stabilization", "Some stabilization must be used in PorousFlowUnsaturated");
80 780 : }
81 :
82 : void
83 750 : PorousFlowUnsaturated::addMaterialDependencies()
84 : {
85 750 : PorousFlowSinglePhaseBase::addMaterialDependencies();
86 :
87 : // Add necessary objects to list of PorousFlow objects added by this action
88 750 : _included_objects.push_back("PorousFlowAdvectiveFlux");
89 :
90 750 : if (_transient)
91 1370 : _included_objects.push_back("PorousFlowMassTimeDerivative");
92 :
93 750 : if (_mechanical && _transient)
94 270 : _included_objects.push_back("PorousFlowMassVolumetricExpansion");
95 :
96 750 : if (_thermal)
97 290 : _included_objects.push_back("PorousFlowHeatAdvection");
98 :
99 750 : if (_add_saturation_aux)
100 1410 : _included_objects.push_back("SaturationAux");
101 :
102 750 : if (_stabilization == StabilizationEnum::KT)
103 90 : _included_objects.push_back("PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent");
104 :
105 750 : if (_stabilization == StabilizationEnum::KT && _thermal)
106 90 : _included_objects.push_back("PorousFlowAdvectiveFluxCalculatorUnsaturatedHeat");
107 750 : }
108 :
109 : void
110 146 : PorousFlowUnsaturated::addKernels()
111 : {
112 146 : PorousFlowSinglePhaseBase::addKernels();
113 :
114 : // add the kernels
115 146 : if (_stabilization == StabilizationEnum::Full)
116 : {
117 137 : const std::string kernel_type = "PorousFlowAdvectiveFlux";
118 137 : InputParameters params = _factory.getValidParams(kernel_type);
119 137 : if (_subdomain_names_set)
120 54 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
121 274 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
122 137 : params.set<RealVectorValue>("gravity") = _gravity;
123 :
124 148 : for (unsigned i = 0; i < _num_mass_fraction_vars; ++i)
125 : {
126 11 : const std::string kernel_name = "PorousFlowUnsaturated_AdvectiveFlux" + Moose::stringify(i);
127 11 : params.set<unsigned int>("fluid_component") = i;
128 22 : params.set<NonlinearVariableName>("variable") = _mass_fraction_vars[i];
129 11 : _problem->addKernel(kernel_type, kernel_name, params);
130 : }
131 : const std::string kernel_name =
132 137 : "PorousFlowUnsaturated_AdvectiveFlux" + Moose::stringify(_num_mass_fraction_vars);
133 137 : params.set<unsigned int>("fluid_component") = _num_mass_fraction_vars;
134 274 : params.set<NonlinearVariableName>("variable") = _pp_var;
135 137 : _problem->addKernel(kernel_type, kernel_name, params);
136 137 : }
137 9 : else if (_stabilization == StabilizationEnum::KT)
138 : {
139 9 : const std::string kernel_type = "PorousFlowFluxLimitedTVDAdvection";
140 9 : InputParameters params = _factory.getValidParams(kernel_type);
141 9 : if (_subdomain_names_set)
142 18 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
143 18 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
144 :
145 9 : for (unsigned i = 0; i < _num_mass_fraction_vars; ++i)
146 : {
147 0 : const std::string kernel_name = "PorousFlowFluxLimited_DarcyFlow" + Moose::stringify(i);
148 0 : params.set<UserObjectName>("advective_flux_calculator") =
149 0 : "PorousFlowUnsaturated_AC_" + Moose::stringify(i);
150 0 : params.set<NonlinearVariableName>("variable") = _mass_fraction_vars[i];
151 0 : _problem->addKernel(kernel_type, kernel_name, params);
152 : }
153 : const std::string kernel_name =
154 9 : "PorousFlowFluxLimited_DarcyFlow" + Moose::stringify(_num_mass_fraction_vars);
155 18 : params.set<NonlinearVariableName>("variable") = _pp_var;
156 18 : params.set<UserObjectName>("advective_flux_calculator") =
157 9 : "PorousFlowUnsaturated_AC_" + Moose::stringify(_num_mass_fraction_vars);
158 9 : _problem->addKernel(kernel_type, kernel_name, params);
159 9 : }
160 :
161 146 : if (_transient)
162 : {
163 137 : std::string kernel_name = "PorousFlowUnsaturated_MassTimeDerivative";
164 137 : std::string kernel_type = "PorousFlowMassTimeDerivative";
165 137 : InputParameters params = _factory.getValidParams(kernel_type);
166 137 : if (_subdomain_names_set)
167 54 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
168 274 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
169 137 : params.set<bool>("strain_at_nearest_qp") = _strain_at_nearest_qp;
170 137 : if (!_base_name.empty())
171 0 : params.set<std::string>("base_name") = _base_name;
172 :
173 148 : for (unsigned i = 0; i < _num_mass_fraction_vars; ++i)
174 : {
175 11 : kernel_name = "PorousFlowUnsaturated_MassTimeDerivative" + Moose::stringify(i);
176 11 : params.set<unsigned int>("fluid_component") = i;
177 22 : params.set<NonlinearVariableName>("variable") = _mass_fraction_vars[i];
178 11 : if (_save_component_rate_in.size() != 0)
179 0 : params.set<std::vector<AuxVariableName>>("save_in") = {_save_component_rate_in[i]};
180 11 : _problem->addKernel(kernel_type, kernel_name, params);
181 : }
182 : kernel_name =
183 137 : "PorousFlowUnsaturated_MassTimeDerivative" + Moose::stringify(_num_mass_fraction_vars);
184 137 : params.set<unsigned int>("fluid_component") = _num_mass_fraction_vars;
185 274 : params.set<NonlinearVariableName>("variable") = _pp_var;
186 137 : if (_save_component_rate_in.size() != 0)
187 18 : params.set<std::vector<AuxVariableName>>("save_in") = {
188 36 : _save_component_rate_in[_num_mass_fraction_vars]};
189 137 : _problem->addKernel(kernel_type, kernel_name, params);
190 137 : }
191 :
192 146 : if (_mechanical && _transient)
193 : {
194 27 : std::string kernel_name = "PorousFlowUnsaturated_MassVolumetricExpansion";
195 27 : std::string kernel_type = "PorousFlowMassVolumetricExpansion";
196 27 : InputParameters params = _factory.getValidParams(kernel_type);
197 27 : if (_subdomain_names_set)
198 54 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
199 54 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
200 27 : params.set<bool>("strain_at_nearest_qp") = _strain_at_nearest_qp;
201 :
202 27 : for (unsigned i = 0; i < _num_mass_fraction_vars; ++i)
203 : {
204 0 : kernel_name = "PorousFlowUnsaturated_MassVolumetricExpansion" + Moose::stringify(i);
205 0 : params.set<unsigned>("fluid_component") = i;
206 0 : params.set<NonlinearVariableName>("variable") = _mass_fraction_vars[i];
207 0 : _problem->addKernel(kernel_type, kernel_name, params);
208 : }
209 : kernel_name =
210 27 : "PorousFlowUnsaturated_MassVolumetricExpansion" + Moose::stringify(_num_mass_fraction_vars);
211 27 : params.set<unsigned>("fluid_component") = _num_mass_fraction_vars;
212 54 : params.set<NonlinearVariableName>("variable") = _pp_var;
213 27 : _problem->addKernel(kernel_type, kernel_name, params);
214 27 : }
215 :
216 146 : if (_thermal)
217 : {
218 29 : if (_stabilization == StabilizationEnum::Full)
219 : {
220 20 : const std::string kernel_name = "PorousFlowUnsaturated_HeatAdvection";
221 20 : const std::string kernel_type = "PorousFlowHeatAdvection";
222 20 : InputParameters params = _factory.getValidParams(kernel_type);
223 20 : if (_subdomain_names_set)
224 18 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
225 40 : params.set<NonlinearVariableName>("variable") = _temperature_var[0];
226 40 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
227 20 : params.set<RealVectorValue>("gravity") = _gravity;
228 20 : _problem->addKernel(kernel_type, kernel_name, params);
229 20 : }
230 9 : else if (_stabilization == StabilizationEnum::KT)
231 : {
232 9 : const std::string kernel_name = "PorousFlowUnsaturated_HeatAdvection";
233 9 : const std::string kernel_type = "PorousFlowFluxLimitedTVDAdvection";
234 9 : InputParameters params = _factory.getValidParams(kernel_type);
235 9 : if (_subdomain_names_set)
236 18 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
237 18 : params.set<NonlinearVariableName>("variable") = _temperature_var[0];
238 18 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
239 18 : params.set<UserObjectName>("advective_flux_calculator") = "PorousFlowUnsaturatedHeat_AC";
240 9 : _problem->addKernel(kernel_type, kernel_name, params);
241 9 : }
242 : }
243 146 : }
244 :
245 : void
246 156 : PorousFlowUnsaturated::addUserObjects()
247 : {
248 156 : PorousFlowSinglePhaseBase::addUserObjects();
249 :
250 : // Add the capillary pressure UserObject
251 156 : addCapillaryPressureVG(_van_genuchten_m, _van_genuchten_alpha, _capillary_pressure_name, _s_res);
252 :
253 : // add Advective Flux calculator UserObjects, if required
254 156 : if (_stabilization == StabilizationEnum::KT)
255 : {
256 9 : for (unsigned i = 0; i < _num_mass_fraction_vars; ++i)
257 : {
258 0 : const std::string userobject_name = "PorousFlowUnsaturated_AC_" + Moose::stringify(i);
259 0 : addAdvectiveFluxCalculatorUnsaturatedMultiComponent(0, i, true, userobject_name);
260 : }
261 : const std::string userobject_name =
262 9 : "PorousFlowUnsaturated_AC_" + Moose::stringify(_num_mass_fraction_vars);
263 9 : if (_num_mass_fraction_vars == 0)
264 18 : addAdvectiveFluxCalculatorUnsaturated(0, true, userobject_name); // 1 component only
265 : else
266 0 : addAdvectiveFluxCalculatorUnsaturatedMultiComponent(
267 0 : 0, _num_mass_fraction_vars, true, userobject_name);
268 :
269 9 : if (_thermal)
270 : {
271 9 : const std::string userobject_name = "PorousFlowUnsaturatedHeat_AC";
272 18 : addAdvectiveFluxCalculatorUnsaturatedHeat(0, true, userobject_name);
273 : }
274 : }
275 156 : }
276 :
277 : void
278 146 : PorousFlowUnsaturated::addMaterials()
279 : {
280 146 : PorousFlowSinglePhaseBase::addMaterials();
281 :
282 292 : if (_deps.dependsOn(_included_objects, "pressure_saturation_qp"))
283 : {
284 146 : const std::string material_type = "PorousFlow1PhaseP";
285 146 : const std::string material_name = "PorousFlowUnsaturated_1PhaseP_VG_qp";
286 146 : InputParameters params = _factory.getValidParams(material_type);
287 146 : if (_subdomain_names_set)
288 72 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
289 292 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
290 438 : params.set<std::vector<VariableName>>("porepressure") = {_pp_var};
291 292 : params.set<UserObjectName>("capillary_pressure") = _capillary_pressure_name;
292 146 : params.set<bool>("at_nodes") = false;
293 146 : _problem->addMaterial(material_type, material_name, params);
294 146 : }
295 292 : if (_deps.dependsOn(_included_objects, "pressure_saturation_nodal"))
296 : {
297 146 : const std::string material_type = "PorousFlow1PhaseP";
298 146 : const std::string material_name = "PorousFlowUnsaturated_1PhaseP_VG_nodal";
299 146 : InputParameters params = _factory.getValidParams(material_type);
300 146 : if (_subdomain_names_set)
301 72 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
302 292 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
303 438 : params.set<std::vector<VariableName>>("porepressure") = {_pp_var};
304 292 : params.set<UserObjectName>("capillary_pressure") = _capillary_pressure_name;
305 146 : params.set<bool>("at_nodes") = true;
306 146 : _problem->addMaterial(material_type, material_name, params);
307 146 : }
308 :
309 292 : if (_deps.dependsOn(_included_objects, "relative_permeability_qp"))
310 : {
311 126 : if (_relperm_type == RelpermTypeChoiceEnum::FLAC)
312 90 : addRelativePermeabilityFLAC(false, 0, _relative_permeability_exponent, _s_res, _s_res);
313 : else
314 36 : addRelativePermeabilityCorey(false, 0, _relative_permeability_exponent, _s_res, _s_res);
315 : }
316 :
317 292 : if (_deps.dependsOn(_included_objects, "relative_permeability_nodal"))
318 : {
319 146 : if (_relperm_type == RelpermTypeChoiceEnum::FLAC)
320 101 : addRelativePermeabilityFLAC(true, 0, _relative_permeability_exponent, _s_res, _s_res);
321 : else
322 45 : addRelativePermeabilityCorey(true, 0, _relative_permeability_exponent, _s_res, _s_res);
323 : }
324 :
325 292 : if (_deps.dependsOn(_included_objects, "volumetric_strain_qp") ||
326 265 : _deps.dependsOn(_included_objects, "volumetric_strain_nodal"))
327 27 : addVolumetricStrainMaterial(_coupled_displacements, _base_name);
328 146 : }
329 :
330 : void
331 302 : PorousFlowUnsaturated::addAuxObjects()
332 : {
333 302 : PorousFlowSinglePhaseBase::addAuxObjects();
334 :
335 302 : if (_add_saturation_aux)
336 284 : addSaturationAux(0);
337 302 : }
|
