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 "PorousFlowOutflowBC.h"
11 :
12 : #include "MooseVariable.h"
13 : #include "libmesh/string_to_enum.h"
14 : #include "libmesh/quadrature.h"
15 :
16 : registerMooseObject("PorousFlowApp", PorousFlowOutflowBC);
17 :
18 : InputParameters
19 383 : PorousFlowOutflowBC::validParams()
20 : {
21 383 : InputParameters params = IntegratedBC::validParams();
22 766 : params.addRequiredParam<UserObjectName>(
23 : "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
24 766 : params.addRequiredParam<RealVectorValue>("gravity",
25 : "Gravitational acceleration vector downwards (m/s^2)");
26 766 : MooseEnum flux_type_enum("fluid heat", "fluid");
27 766 : params.addParam<MooseEnum>(
28 : "flux_type",
29 : flux_type_enum,
30 : "The type of boundary condition to apply. 'fluid' means this boundary condition will allow "
31 : "a fluid component to flow freely from the boundary. 'heat' means this boundary condition "
32 : "will allow heat energy to flow freely from the boundary");
33 766 : params.addParam<unsigned int>(
34 : "mass_fraction_component",
35 766 : 0,
36 : "The index corresponding to a fluid "
37 : "component. If supplied, the residual contribution will be "
38 : "multiplied by the mass fraction, corresponding to allowing the "
39 : "given mass fraction to flow freely from the boundary. This is ignored if flux_type = heat");
40 766 : params.addParam<bool>("multiply_by_density",
41 766 : true,
42 : "If true, this BC represents mass flux. If false, it represents volume "
43 : "flux. User input of this flag is ignored if flux_type = heat as "
44 : "multiply_by_density should always be true in that case");
45 766 : params.addParam<bool>(
46 : "include_relperm",
47 766 : true,
48 : "If true, the Darcy flux will include the relative permeability. If false, the relative "
49 : "permeability will not be used, which must only be used for fully-saturated situations "
50 : "where there is no notion of relative permeability");
51 766 : params.addParam<Real>(
52 : "multiplier",
53 766 : 1.0,
54 : "Multiply the flux by this number. This is mainly used for testing purposes");
55 766 : params.addParamNamesToGroup("multiplier", "Advanced");
56 383 : params.addClassDescription(
57 : "Applies an 'outflow' boundary condition, which allows fluid components or heat energy to "
58 : "flow freely out of the boundary as if it weren't there. This is fully upwinded");
59 383 : return params;
60 383 : }
61 :
62 216 : PorousFlowOutflowBC::PorousFlowOutflowBC(const InputParameters & parameters)
63 : : IntegratedBC(parameters),
64 216 : _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
65 216 : _num_phases(_dictator.numPhases()),
66 216 : _perm_derivs(_dictator.usePermDerivs()),
67 432 : _flux_type(getParam<MooseEnum>("flux_type").getEnum<FluxTypeChoiceEnum>()),
68 432 : _sp(getParam<unsigned>("mass_fraction_component")),
69 216 : _multiply_by_density(
70 398 : _flux_type == FluxTypeChoiceEnum::FLUID ? getParam<bool>("multiply_by_density") : true),
71 432 : _include_relperm(getParam<bool>("include_relperm")),
72 432 : _gravity(getParam<RealVectorValue>("gravity")),
73 432 : _multiplier(getParam<Real>("multiplier")),
74 432 : _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
75 432 : _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
76 : "dPorousFlow_grad_porepressure_qp_dgradvar")),
77 432 : _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
78 : "dPorousFlow_grad_porepressure_qp_dvar")),
79 432 : _fluid_density_qp(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
80 432 : _dfluid_density_qp_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
81 : "dPorousFlow_fluid_phase_density_qp_dvar")),
82 432 : _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
83 216 : _dpermeability_dvar(
84 216 : getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
85 432 : _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
86 : "dPorousFlow_permeability_qp_dgradvar")),
87 432 : _fluid_viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")),
88 216 : _dfluid_viscosity_dvar(
89 216 : getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
90 :
91 216 : _has_density(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
92 430 : hasMaterialProperty<std::vector<std::vector<Real>>>(
93 : "dPorousFlow_fluid_phase_density_nodal_dvar")),
94 216 : _has_mass_fraction(
95 216 : hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
96 430 : hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
97 : "dPorousFlow_mass_frac_nodal_dvar")),
98 216 : _has_relperm(hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
99 430 : hasMaterialProperty<std::vector<std::vector<Real>>>(
100 : "dPorousFlow_relative_permeability_nodal_dvar")),
101 216 : _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
102 428 : hasMaterialProperty<std::vector<std::vector<Real>>>(
103 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
104 216 : _has_thermal_conductivity(
105 216 : hasMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp") &&
106 263 : hasMaterialProperty<std::vector<RealTensorValue>>(
107 : "dPorousFlow_thermal_conductivity_qp_dvar")),
108 430 : _has_t(hasMaterialProperty<RealGradient>("PorousFlow_grad_temperature_qp") &&
109 432 : hasMaterialProperty<std::vector<RealGradient>>("dPorousFlow_grad_temperature_qp_dvar") &&
110 430 : hasMaterialProperty<std::vector<Real>>("dPorousFlow_grad_temperature_qp_dgradvar")),
111 :
112 430 : _fluid_density_node(_has_density ? &getMaterialProperty<std::vector<Real>>(
113 : "PorousFlow_fluid_phase_density_nodal")
114 : : nullptr),
115 430 : _dfluid_density_node_dvar(_has_density ? &getMaterialProperty<std::vector<std::vector<Real>>>(
116 : "dPorousFlow_fluid_phase_density_nodal_dvar")
117 : : nullptr),
118 430 : _relative_permeability(_has_relperm ? &getMaterialProperty<std::vector<Real>>(
119 : "PorousFlow_relative_permeability_nodal")
120 : : nullptr),
121 432 : _drelative_permeability_dvar(_has_relperm
122 430 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
123 : "dPorousFlow_relative_permeability_nodal_dvar")
124 : : nullptr),
125 430 : _mass_fractions(_has_mass_fraction ? &getMaterialProperty<std::vector<std::vector<Real>>>(
126 : "PorousFlow_mass_frac_nodal")
127 : : nullptr),
128 432 : _dmass_fractions_dvar(_has_mass_fraction
129 430 : ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
130 : "dPorousFlow_mass_frac_nodal_dvar")
131 : : nullptr),
132 428 : _enthalpy(_has_enthalpy ? &getMaterialPropertyByName<std::vector<Real>>(
133 : "PorousFlow_fluid_phase_enthalpy_nodal")
134 : : nullptr),
135 428 : _denthalpy_dvar(_has_enthalpy ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
136 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
137 : : nullptr),
138 263 : _thermal_conductivity(_has_thermal_conductivity ? &getMaterialProperty<RealTensorValue>(
139 : "PorousFlow_thermal_conductivity_qp")
140 : : nullptr),
141 432 : _dthermal_conductivity_dvar(_has_thermal_conductivity
142 263 : ? &getMaterialProperty<std::vector<RealTensorValue>>(
143 : "dPorousFlow_thermal_conductivity_qp_dvar")
144 : : nullptr),
145 430 : _grad_t(_has_t ? &getMaterialProperty<RealGradient>("PorousFlow_grad_temperature_qp")
146 : : nullptr),
147 430 : _dgrad_t_dvar(_has_t ? &getMaterialProperty<std::vector<RealGradient>>(
148 : "dPorousFlow_grad_temperature_qp_dvar")
149 : : nullptr),
150 216 : _dgrad_t_dgradvar(
151 430 : _has_t ? &getMaterialProperty<std::vector<Real>>("dPorousFlow_grad_temperature_qp_dgradvar")
152 216 : : nullptr)
153 : {
154 216 : if (_flux_type == FluxTypeChoiceEnum::FLUID && _sp >= _dictator.numComponents())
155 2 : paramError("mass_fraction_component",
156 : "The Dictator declares that the maximum fluid component index is ",
157 2 : _dictator.numComponents() - 1,
158 : ", but you have set mass_fraction_component to ",
159 2 : _sp,
160 : ". Remember that indexing starts at 0. Please be assured that the Dictator has "
161 : "noted your error.");
162 :
163 214 : if (_multiply_by_density && !_has_density)
164 2 : mooseError("PorousFlowOutflowBC: You requested to multiply_by_density, but you have no nodal "
165 : "fluid density Material");
166 212 : if (_include_relperm && !_has_relperm)
167 2 : mooseError("PorousFlowOutflowBC: You requested to include the relative permeability, but you "
168 : "have no nodal relative permeability Material");
169 210 : if (_flux_type == FluxTypeChoiceEnum::FLUID && !_has_mass_fraction)
170 2 : mooseError(
171 : "PorousFlowOutflowBC: For flux_type = fluid, you need a nodal mass-fraction Material");
172 208 : if (_flux_type == FluxTypeChoiceEnum::HEAT && !_has_enthalpy)
173 2 : mooseError(
174 : "PorousFlowOutflowBC: For flux_type = heat, you need a nodal fluid enthalpy Material");
175 206 : if (_flux_type == FluxTypeChoiceEnum::HEAT && !_has_thermal_conductivity)
176 2 : mooseError(
177 : "PorousFlowOutflowBC: For flux_type = heat, you need a thermal conductivity Material");
178 204 : if (_flux_type == FluxTypeChoiceEnum::HEAT && !_has_t)
179 2 : mooseError("PorousFlowOutflowBC: For flux_type = heat, you need a temperature Material");
180 202 : }
181 :
182 : Real
183 854560 : PorousFlowOutflowBC::computeQpResidual()
184 : {
185 : Real advective_term = 0.0;
186 1893560 : for (unsigned ph = 0; ph < _num_phases; ++ph)
187 1039000 : advective_term -= darcy(ph) * mobility(ph) * prefactor(ph);
188 :
189 854560 : if (_flux_type == FluxTypeChoiceEnum::FLUID)
190 615536 : return _test[_i][_qp] * advective_term * _multiplier;
191 :
192 239024 : const Real conduction_term = -_normals[_qp] * ((*_thermal_conductivity)[_qp] * (*_grad_t)[_qp]);
193 239024 : return _test[_i][_qp] * (conduction_term + advective_term) * _multiplier;
194 : }
195 :
196 : Real
197 399452 : PorousFlowOutflowBC::computeQpJacobian()
198 : {
199 399452 : return jac(_var.number());
200 : }
201 :
202 : Real
203 414940 : PorousFlowOutflowBC::computeQpOffDiagJacobian(unsigned int jvar)
204 : {
205 414940 : return jac(jvar);
206 : }
207 :
208 : Real
209 814392 : PorousFlowOutflowBC::jac(unsigned int jvar) const
210 : {
211 814392 : if (_dictator.notPorousFlowVariable(jvar))
212 : return 0.0;
213 814392 : const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
214 :
215 : // For _i != _j, note:
216 : // since the only non-upwinded contribution to the residual is
217 : // from Darcy and thermal_conductivity, the only contribution
218 : // of the residual at node _i from changing jvar at node _j is through
219 : // the derivative of Darcy or thermal_conductivity
220 :
221 : Real advective_term_prime = 0.0;
222 1776560 : for (unsigned ph = 0; ph < _num_phases; ++ph)
223 : {
224 : Real darcyprime =
225 962168 : _normals[_qp] *
226 962168 : (_permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] +
227 962168 : _dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp] -
228 962168 : _dfluid_density_qp_dvar[_qp][ph][pvar] * _phi[_j][_qp] * _gravity));
229 962168 : if (_perm_derivs)
230 : {
231 0 : darcyprime += _normals[_qp] * (_dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
232 0 : (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
233 0 : for (const auto i : make_range(Moose::dim))
234 0 : darcyprime +=
235 0 : _normals[_qp] * (_dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
236 0 : (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
237 : }
238 962168 : if (_i != _j)
239 820892 : advective_term_prime -= prefactor(ph) * mobility(ph) * darcyprime;
240 : else
241 : {
242 141276 : const Real dar = darcy(ph);
243 :
244 141276 : Real mob = 1.0 / _fluid_viscosity[_i][ph];
245 : Real mobprime =
246 141276 : -_dfluid_viscosity_dvar[_i][ph][pvar] / Utility::pow<2>(_fluid_viscosity[_i][ph]);
247 141276 : if (_multiply_by_density)
248 : {
249 141276 : const Real densityprime = (*_dfluid_density_node_dvar)[_i][ph][pvar];
250 141276 : mobprime = densityprime * mob + (*_fluid_density_node)[_i][ph] * mobprime;
251 141276 : mob *= (*_fluid_density_node)[_i][ph];
252 : }
253 141276 : if (_include_relperm)
254 : {
255 137376 : const Real relpermprime = (*_drelative_permeability_dvar)[_i][ph][pvar];
256 137376 : mobprime = relpermprime * mob + (*_relative_permeability)[_i][ph] * mobprime;
257 137376 : mob *= (*_relative_permeability)[_i][ph];
258 : }
259 :
260 141276 : const Real pre = prefactor(ph);
261 : const Real preprime =
262 141276 : (_flux_type == FluxTypeChoiceEnum::FLUID ? (*_dmass_fractions_dvar)[_i][ph][_sp][pvar]
263 44160 : : (*_denthalpy_dvar)[_i][ph][pvar]);
264 :
265 141276 : advective_term_prime -= preprime * mob * dar + pre * mobprime * dar + pre * mob * darcyprime;
266 : }
267 : }
268 814392 : if (_flux_type == FluxTypeChoiceEnum::FLUID)
269 545592 : return _test[_i][_qp] * advective_term_prime * _multiplier;
270 :
271 : const Real conduction_term_prime =
272 268800 : -_normals[_qp] *
273 268800 : ((*_dthermal_conductivity_dvar)[_qp][pvar] * (*_grad_t)[_qp] +
274 268800 : (*_thermal_conductivity)[_qp] * (*_dgrad_t_dvar)[_qp][pvar]) *
275 268800 : _phi[_j][_qp] -
276 : _normals[_qp] *
277 268800 : ((*_thermal_conductivity)[_qp] * (*_dgrad_t_dgradvar)[_qp][pvar] * _grad_phi[_j][_qp]);
278 268800 : return _test[_i][_qp] * (conduction_term_prime + advective_term_prime) * _multiplier;
279 : }
280 :
281 : Real
282 1180276 : PorousFlowOutflowBC::darcy(unsigned ph) const
283 : {
284 1180276 : return _normals[_qp] *
285 1180276 : (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
286 : }
287 :
288 : Real
289 1859892 : PorousFlowOutflowBC::mobility(unsigned ph) const
290 : {
291 1859892 : return (_multiply_by_density ? (*_fluid_density_node)[_i][ph] : 1.0) *
292 1859892 : (_include_relperm ? (*_relative_permeability)[_i][ph] : 1.0) / _fluid_viscosity[_i][ph];
293 : }
294 :
295 : Real
296 2001168 : PorousFlowOutflowBC::prefactor(unsigned ph) const
297 : {
298 2001168 : return (_flux_type == FluxTypeChoiceEnum::FLUID ? (*_mass_fractions)[_i][ph][_sp]
299 507824 : : (*_enthalpy)[_i][ph]);
300 : }
|
