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 1584 : PorousFlowUnsaturated::validParams()
28 : {
29 1584 : InputParameters params = PorousFlowSinglePhaseBase::validParams();
30 3168 : params.addParam<bool>("add_saturation_aux", true, "Add an AuxVariable that records saturation");
31 4752 : params.addRangeCheckedParam<Real>(
32 : "van_genuchten_alpha",
33 3168 : 1.0E-6,
34 : "van_genuchten_alpha > 0.0",
35 : "Van Genuchten alpha parameter used to determine saturation from porepressure");
36 4752 : params.addRangeCheckedParam<Real>(
37 : "van_genuchten_m",
38 3168 : 0.6,
39 : "van_genuchten_m > 0 & van_genuchten_m < 1",
40 : "Van Genuchten m parameter used to determine saturation from porepressure");
41 3168 : MooseEnum relperm_type_choice("FLAC Corey", "FLAC");
42 3168 : 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 4752 : params.addRangeCheckedParam<Real>("relative_permeability_exponent",
48 3168 : 3.0,
49 : "relative_permeability_exponent>=0",
50 : "Relative permeability exponent");
51 4752 : params.addRangeCheckedParam<Real>(
52 : "residual_saturation",
53 3168 : 0.0,
54 : "residual_saturation>=0.0 & residual_saturation<1.0",
55 : "Residual saturation to use in the relative permeability expression");
56 1584 : 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 1584 : return params;
65 1584 : }
66 :
67 1584 : PorousFlowUnsaturated::PorousFlowUnsaturated(const InputParameters & params)
68 : : PorousFlowSinglePhaseBase(params),
69 1582 : _add_saturation_aux(getParam<bool>("add_saturation_aux")),
70 3164 : _van_genuchten_alpha(getParam<Real>("van_genuchten_alpha")),
71 3164 : _van_genuchten_m(getParam<Real>("van_genuchten_m")),
72 1582 : _relperm_type(
73 1582 : getParam<MooseEnum>("relative_permeability_type").getEnum<RelpermTypeChoiceEnum>()),
74 3164 : _relative_permeability_exponent(getParam<Real>("relative_permeability_exponent")),
75 3164 : _s_res(getParam<Real>("residual_saturation")),
76 3166 : _capillary_pressure_name("PorousFlowUnsaturated_CapillaryPressureVG")
77 : {
78 1582 : if (_stabilization == StabilizationEnum::None)
79 2 : paramError("stabilization", "Some stabilization must be used in PorousFlowUnsaturated");
80 1580 : }
81 :
82 : void
83 1550 : PorousFlowUnsaturated::addMaterialDependencies()
84 : {
85 1550 : PorousFlowSinglePhaseBase::addMaterialDependencies();
86 :
87 : // Add necessary objects to list of PorousFlow objects added by this action
88 1550 : _included_objects.push_back("PorousFlowAdvectiveFlux");
89 :
90 1550 : if (_transient)
91 2870 : _included_objects.push_back("PorousFlowMassTimeDerivative");
92 :
93 1550 : if (_mechanical && _transient)
94 570 : _included_objects.push_back("PorousFlowMassVolumetricExpansion");
95 :
96 1550 : if (_thermal)
97 590 : _included_objects.push_back("PorousFlowHeatAdvection");
98 :
99 1550 : if (_add_saturation_aux)
100 2910 : _included_objects.push_back("SaturationAux");
101 :
102 1550 : if (_stabilization == StabilizationEnum::KT)
103 190 : _included_objects.push_back("PorousFlowAdvectiveFluxCalculatorUnsaturatedMultiComponent");
104 :
105 1550 : if (_stabilization == StabilizationEnum::KT && _thermal)
106 190 : _included_objects.push_back("PorousFlowAdvectiveFluxCalculatorUnsaturatedHeat");
107 1550 : }
108 :
109 : void
110 306 : PorousFlowUnsaturated::addKernels()
111 : {
112 306 : PorousFlowSinglePhaseBase::addKernels();
113 :
114 : // add the kernels
115 306 : if (_stabilization == StabilizationEnum::Full)
116 : {
117 287 : const std::string kernel_type = "PorousFlowAdvectiveFlux";
118 287 : InputParameters params = _factory.getValidParams(kernel_type);
119 287 : if (_subdomain_names_set)
120 114 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
121 574 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
122 287 : params.set<RealVectorValue>("gravity") = _gravity;
123 :
124 308 : for (unsigned i = 0; i < _num_mass_fraction_vars; ++i)
125 : {
126 21 : const std::string kernel_name = "PorousFlowUnsaturated_AdvectiveFlux" + Moose::stringify(i);
127 21 : params.set<unsigned int>("fluid_component") = i;
128 42 : params.set<NonlinearVariableName>("variable") = _mass_fraction_vars[i];
129 21 : _problem->addKernel(kernel_type, kernel_name, params);
130 : }
131 : const std::string kernel_name =
132 287 : "PorousFlowUnsaturated_AdvectiveFlux" + Moose::stringify(_num_mass_fraction_vars);
133 287 : params.set<unsigned int>("fluid_component") = _num_mass_fraction_vars;
134 574 : params.set<NonlinearVariableName>("variable") = _pp_var;
135 287 : _problem->addKernel(kernel_type, kernel_name, params);
136 287 : }
137 19 : else if (_stabilization == StabilizationEnum::KT)
138 : {
139 19 : const std::string kernel_type = "PorousFlowFluxLimitedTVDAdvection";
140 19 : InputParameters params = _factory.getValidParams(kernel_type);
141 19 : if (_subdomain_names_set)
142 38 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
143 38 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
144 :
145 19 : 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 19 : "PorousFlowFluxLimited_DarcyFlow" + Moose::stringify(_num_mass_fraction_vars);
155 38 : params.set<NonlinearVariableName>("variable") = _pp_var;
156 38 : params.set<UserObjectName>("advective_flux_calculator") =
157 19 : "PorousFlowUnsaturated_AC_" + Moose::stringify(_num_mass_fraction_vars);
158 19 : _problem->addKernel(kernel_type, kernel_name, params);
159 19 : }
160 :
161 306 : if (_transient)
162 : {
163 287 : std::string kernel_name = "PorousFlowUnsaturated_MassTimeDerivative";
164 287 : std::string kernel_type = "PorousFlowMassTimeDerivative";
165 287 : InputParameters params = _factory.getValidParams(kernel_type);
166 287 : if (_subdomain_names_set)
167 114 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
168 574 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
169 287 : params.set<bool>("strain_at_nearest_qp") = _strain_at_nearest_qp;
170 287 : if (!_base_name.empty())
171 0 : params.set<std::string>("base_name") = _base_name;
172 :
173 308 : for (unsigned i = 0; i < _num_mass_fraction_vars; ++i)
174 : {
175 21 : kernel_name = "PorousFlowUnsaturated_MassTimeDerivative" + Moose::stringify(i);
176 21 : params.set<unsigned int>("fluid_component") = i;
177 42 : params.set<NonlinearVariableName>("variable") = _mass_fraction_vars[i];
178 21 : if (_save_component_rate_in.size() != 0)
179 0 : params.set<std::vector<AuxVariableName>>("save_in") = {_save_component_rate_in[i]};
180 21 : _problem->addKernel(kernel_type, kernel_name, params);
181 : }
182 : kernel_name =
183 287 : "PorousFlowUnsaturated_MassTimeDerivative" + Moose::stringify(_num_mass_fraction_vars);
184 287 : params.set<unsigned int>("fluid_component") = _num_mass_fraction_vars;
185 574 : params.set<NonlinearVariableName>("variable") = _pp_var;
186 287 : if (_save_component_rate_in.size() != 0)
187 38 : params.set<std::vector<AuxVariableName>>("save_in") = {
188 76 : _save_component_rate_in[_num_mass_fraction_vars]};
189 287 : _problem->addKernel(kernel_type, kernel_name, params);
190 287 : }
191 :
192 306 : if (_mechanical && _transient)
193 : {
194 57 : std::string kernel_name = "PorousFlowUnsaturated_MassVolumetricExpansion";
195 57 : std::string kernel_type = "PorousFlowMassVolumetricExpansion";
196 57 : InputParameters params = _factory.getValidParams(kernel_type);
197 57 : if (_subdomain_names_set)
198 114 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
199 114 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
200 57 : params.set<bool>("strain_at_nearest_qp") = _strain_at_nearest_qp;
201 :
202 57 : 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 57 : "PorousFlowUnsaturated_MassVolumetricExpansion" + Moose::stringify(_num_mass_fraction_vars);
211 57 : params.set<unsigned>("fluid_component") = _num_mass_fraction_vars;
212 114 : params.set<NonlinearVariableName>("variable") = _pp_var;
213 57 : _problem->addKernel(kernel_type, kernel_name, params);
214 57 : }
215 :
216 306 : if (_thermal)
217 : {
218 59 : if (_stabilization == StabilizationEnum::Full)
219 : {
220 40 : const std::string kernel_name = "PorousFlowUnsaturated_HeatAdvection";
221 40 : const std::string kernel_type = "PorousFlowHeatAdvection";
222 40 : InputParameters params = _factory.getValidParams(kernel_type);
223 40 : if (_subdomain_names_set)
224 38 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
225 80 : params.set<NonlinearVariableName>("variable") = _temperature_var[0];
226 80 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
227 40 : params.set<RealVectorValue>("gravity") = _gravity;
228 40 : _problem->addKernel(kernel_type, kernel_name, params);
229 40 : }
230 19 : else if (_stabilization == StabilizationEnum::KT)
231 : {
232 19 : const std::string kernel_name = "PorousFlowUnsaturated_HeatAdvection";
233 19 : const std::string kernel_type = "PorousFlowFluxLimitedTVDAdvection";
234 19 : InputParameters params = _factory.getValidParams(kernel_type);
235 19 : if (_subdomain_names_set)
236 38 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
237 38 : params.set<NonlinearVariableName>("variable") = _temperature_var[0];
238 38 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
239 38 : params.set<UserObjectName>("advective_flux_calculator") = "PorousFlowUnsaturatedHeat_AC";
240 19 : _problem->addKernel(kernel_type, kernel_name, params);
241 19 : }
242 : }
243 306 : }
244 :
245 : void
246 316 : PorousFlowUnsaturated::addUserObjects()
247 : {
248 316 : PorousFlowSinglePhaseBase::addUserObjects();
249 :
250 : // Add the capillary pressure UserObject
251 316 : addCapillaryPressureVG(_van_genuchten_m, _van_genuchten_alpha, _capillary_pressure_name);
252 :
253 : // add Advective Flux calculator UserObjects, if required
254 316 : if (_stabilization == StabilizationEnum::KT)
255 : {
256 19 : 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 19 : "PorousFlowUnsaturated_AC_" + Moose::stringify(_num_mass_fraction_vars);
263 19 : if (_num_mass_fraction_vars == 0)
264 38 : 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 19 : if (_thermal)
270 : {
271 19 : const std::string userobject_name = "PorousFlowUnsaturatedHeat_AC";
272 38 : addAdvectiveFluxCalculatorUnsaturatedHeat(0, true, userobject_name);
273 : }
274 : }
275 316 : }
276 :
277 : void
278 306 : PorousFlowUnsaturated::addMaterials()
279 : {
280 306 : PorousFlowSinglePhaseBase::addMaterials();
281 :
282 612 : if (_deps.dependsOn(_included_objects, "pressure_saturation_qp"))
283 : {
284 306 : const std::string material_type = "PorousFlow1PhaseP";
285 306 : const std::string material_name = "PorousFlowUnsaturated_1PhaseP_VG_qp";
286 306 : InputParameters params = _factory.getValidParams(material_type);
287 306 : if (_subdomain_names_set)
288 152 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
289 612 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
290 918 : params.set<std::vector<VariableName>>("porepressure") = {_pp_var};
291 612 : params.set<UserObjectName>("capillary_pressure") = _capillary_pressure_name;
292 306 : params.set<bool>("at_nodes") = false;
293 306 : _problem->addMaterial(material_type, material_name, params);
294 306 : }
295 612 : if (_deps.dependsOn(_included_objects, "pressure_saturation_nodal"))
296 : {
297 306 : const std::string material_type = "PorousFlow1PhaseP";
298 306 : const std::string material_name = "PorousFlowUnsaturated_1PhaseP_VG_nodal";
299 306 : InputParameters params = _factory.getValidParams(material_type);
300 306 : if (_subdomain_names_set)
301 152 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
302 612 : params.set<UserObjectName>("PorousFlowDictator") = _dictator_name;
303 918 : params.set<std::vector<VariableName>>("porepressure") = {_pp_var};
304 612 : params.set<UserObjectName>("capillary_pressure") = _capillary_pressure_name;
305 306 : params.set<bool>("at_nodes") = true;
306 306 : _problem->addMaterial(material_type, material_name, params);
307 306 : }
308 :
309 612 : if (_deps.dependsOn(_included_objects, "relative_permeability_qp"))
310 : {
311 285 : if (_relperm_type == RelpermTypeChoiceEnum::FLAC)
312 209 : addRelativePermeabilityFLAC(false, 0, _relative_permeability_exponent, _s_res, _s_res);
313 : else
314 76 : addRelativePermeabilityCorey(false, 0, _relative_permeability_exponent, _s_res, _s_res);
315 : }
316 :
317 612 : if (_deps.dependsOn(_included_objects, "relative_permeability_nodal"))
318 : {
319 306 : if (_relperm_type == RelpermTypeChoiceEnum::FLAC)
320 211 : addRelativePermeabilityFLAC(true, 0, _relative_permeability_exponent, _s_res, _s_res);
321 : else
322 95 : addRelativePermeabilityCorey(true, 0, _relative_permeability_exponent, _s_res, _s_res);
323 : }
324 :
325 612 : if (_deps.dependsOn(_included_objects, "volumetric_strain_qp") ||
326 555 : _deps.dependsOn(_included_objects, "volumetric_strain_nodal"))
327 57 : addVolumetricStrainMaterial(_coupled_displacements, _base_name);
328 306 : }
329 :
330 : void
331 622 : PorousFlowUnsaturated::addAuxObjects()
332 : {
333 622 : PorousFlowSinglePhaseBase::addAuxObjects();
334 :
335 622 : if (_add_saturation_aux)
336 584 : addSaturationAux(0);
337 622 : }
|
