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 217 : PorousFlowOutflowBC::validParams()
20 : {
21 217 : InputParameters params = IntegratedBC::validParams();
22 434 : params.addRequiredParam<UserObjectName>(
23 : "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
24 434 : params.addRequiredParam<RealVectorValue>("gravity",
25 : "Gravitational acceleration vector downwards (m/s^2)");
26 434 : MooseEnum flux_type_enum("fluid heat", "fluid");
27 434 : 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 434 : params.addParam<unsigned int>(
34 : "mass_fraction_component",
35 434 : 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 434 : params.addParam<bool>("multiply_by_density",
41 434 : 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 434 : params.addParam<bool>(
46 : "include_relperm",
47 434 : 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 434 : params.addParam<Real>(
52 : "multiplier",
53 434 : 1.0,
54 : "Multiply the flux by this number. This is mainly used for testing purposes");
55 434 : params.addParamNamesToGroup("multiplier", "Advanced");
56 217 : 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 217 : return params;
60 217 : }
61 :
62 124 : PorousFlowOutflowBC::PorousFlowOutflowBC(const InputParameters & parameters)
63 : : IntegratedBC(parameters),
64 124 : _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
65 124 : _num_phases(_dictator.numPhases()),
66 124 : _perm_derivs(_dictator.usePermDerivs()),
67 248 : _flux_type(getParam<MooseEnum>("flux_type").getEnum<FluxTypeChoiceEnum>()),
68 248 : _sp(getParam<unsigned>("mass_fraction_component")),
69 124 : _multiply_by_density(
70 226 : _flux_type == FluxTypeChoiceEnum::FLUID ? getParam<bool>("multiply_by_density") : true),
71 248 : _include_relperm(getParam<bool>("include_relperm")),
72 248 : _gravity(getParam<RealVectorValue>("gravity")),
73 248 : _multiplier(getParam<Real>("multiplier")),
74 248 : _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
75 248 : _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
76 : "dPorousFlow_grad_porepressure_qp_dgradvar")),
77 248 : _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
78 : "dPorousFlow_grad_porepressure_qp_dvar")),
79 248 : _fluid_density_qp(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
80 248 : _dfluid_density_qp_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
81 : "dPorousFlow_fluid_phase_density_qp_dvar")),
82 248 : _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
83 124 : _dpermeability_dvar(
84 124 : getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
85 248 : _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
86 : "dPorousFlow_permeability_qp_dgradvar")),
87 248 : _fluid_viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")),
88 124 : _dfluid_viscosity_dvar(
89 124 : getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
90 :
91 124 : _has_density(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
92 246 : hasMaterialProperty<std::vector<std::vector<Real>>>(
93 : "dPorousFlow_fluid_phase_density_nodal_dvar")),
94 124 : _has_mass_fraction(
95 124 : hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
96 246 : hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
97 : "dPorousFlow_mass_frac_nodal_dvar")),
98 124 : _has_relperm(hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
99 246 : hasMaterialProperty<std::vector<std::vector<Real>>>(
100 : "dPorousFlow_relative_permeability_nodal_dvar")),
101 124 : _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
102 244 : hasMaterialProperty<std::vector<std::vector<Real>>>(
103 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
104 124 : _has_thermal_conductivity(
105 124 : hasMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp") &&
106 159 : hasMaterialProperty<std::vector<RealTensorValue>>(
107 : "dPorousFlow_thermal_conductivity_qp_dvar")),
108 246 : _has_t(hasMaterialProperty<RealGradient>("PorousFlow_grad_temperature_qp") &&
109 248 : hasMaterialProperty<std::vector<RealGradient>>("dPorousFlow_grad_temperature_qp_dvar") &&
110 246 : hasMaterialProperty<std::vector<Real>>("dPorousFlow_grad_temperature_qp_dgradvar")),
111 :
112 246 : _fluid_density_node(_has_density ? &getMaterialProperty<std::vector<Real>>(
113 : "PorousFlow_fluid_phase_density_nodal")
114 : : nullptr),
115 246 : _dfluid_density_node_dvar(_has_density ? &getMaterialProperty<std::vector<std::vector<Real>>>(
116 : "dPorousFlow_fluid_phase_density_nodal_dvar")
117 : : nullptr),
118 246 : _relative_permeability(_has_relperm ? &getMaterialProperty<std::vector<Real>>(
119 : "PorousFlow_relative_permeability_nodal")
120 : : nullptr),
121 248 : _drelative_permeability_dvar(_has_relperm
122 246 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
123 : "dPorousFlow_relative_permeability_nodal_dvar")
124 : : nullptr),
125 246 : _mass_fractions(_has_mass_fraction ? &getMaterialProperty<std::vector<std::vector<Real>>>(
126 : "PorousFlow_mass_frac_nodal")
127 : : nullptr),
128 248 : _dmass_fractions_dvar(_has_mass_fraction
129 246 : ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
130 : "dPorousFlow_mass_frac_nodal_dvar")
131 : : nullptr),
132 244 : _enthalpy(_has_enthalpy ? &getMaterialPropertyByName<std::vector<Real>>(
133 : "PorousFlow_fluid_phase_enthalpy_nodal")
134 : : nullptr),
135 244 : _denthalpy_dvar(_has_enthalpy ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
136 : "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
137 : : nullptr),
138 159 : _thermal_conductivity(_has_thermal_conductivity ? &getMaterialProperty<RealTensorValue>(
139 : "PorousFlow_thermal_conductivity_qp")
140 : : nullptr),
141 248 : _dthermal_conductivity_dvar(_has_thermal_conductivity
142 159 : ? &getMaterialProperty<std::vector<RealTensorValue>>(
143 : "dPorousFlow_thermal_conductivity_qp_dvar")
144 : : nullptr),
145 246 : _grad_t(_has_t ? &getMaterialProperty<RealGradient>("PorousFlow_grad_temperature_qp")
146 : : nullptr),
147 246 : _dgrad_t_dvar(_has_t ? &getMaterialProperty<std::vector<RealGradient>>(
148 : "dPorousFlow_grad_temperature_qp_dvar")
149 : : nullptr),
150 124 : _dgrad_t_dgradvar(
151 246 : _has_t ? &getMaterialProperty<std::vector<Real>>("dPorousFlow_grad_temperature_qp_dgradvar")
152 124 : : nullptr)
153 : {
154 124 : 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 122 : 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 120 : 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 118 : 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 116 : if (_flux_type == FluxTypeChoiceEnum::HEAT && !_has_enthalpy)
173 2 : mooseError(
174 : "PorousFlowOutflowBC: For flux_type = heat, you need a nodal fluid enthalpy Material");
175 114 : if (_flux_type == FluxTypeChoiceEnum::HEAT && !_has_thermal_conductivity)
176 2 : mooseError(
177 : "PorousFlowOutflowBC: For flux_type = heat, you need a thermal conductivity Material");
178 112 : if (_flux_type == FluxTypeChoiceEnum::HEAT && !_has_t)
179 2 : mooseError("PorousFlowOutflowBC: For flux_type = heat, you need a temperature Material");
180 110 : }
181 :
182 : Real
183 840360 : PorousFlowOutflowBC::computeQpResidual()
184 : {
185 : Real advective_term = 0.0;
186 1865124 : for (unsigned ph = 0; ph < _num_phases; ++ph)
187 1024764 : advective_term -= darcy(ph) * mobility(ph) * prefactor(ph);
188 :
189 840360 : if (_flux_type == FluxTypeChoiceEnum::FLUID)
190 607016 : return _test[_i][_qp] * advective_term * _multiplier;
191 :
192 233344 : const Real conduction_term = -_normals[_qp] * ((*_thermal_conductivity)[_qp] * (*_grad_t)[_qp]);
193 233344 : return _test[_i][_qp] * (conduction_term + advective_term) * _multiplier;
194 : }
195 :
196 : Real
197 365716 : PorousFlowOutflowBC::computeQpJacobian()
198 : {
199 365716 : return jac(_var.number());
200 : }
201 :
202 : Real
203 400404 : PorousFlowOutflowBC::computeQpOffDiagJacobian(unsigned int jvar)
204 : {
205 400404 : return jac(jvar);
206 : }
207 :
208 : Real
209 766120 : PorousFlowOutflowBC::jac(unsigned int jvar) const
210 : {
211 766120 : if (_dictator.notPorousFlowVariable(jvar))
212 : return 0.0;
213 766120 : 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 1679920 : for (unsigned ph = 0; ph < _num_phases; ++ph)
223 : {
224 : Real darcyprime =
225 913800 : _normals[_qp] *
226 913800 : (_permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] +
227 913800 : _dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp] -
228 913800 : _dfluid_density_qp_dvar[_qp][ph][pvar] * _phi[_j][_qp] * _gravity));
229 913800 : 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 913800 : if (_i != _j)
239 785348 : advective_term_prime -= prefactor(ph) * mobility(ph) * darcyprime;
240 : else
241 : {
242 128452 : const Real dar = darcy(ph);
243 :
244 128452 : Real mob = 1.0 / _fluid_viscosity[_i][ph];
245 : Real mobprime =
246 128452 : -_dfluid_viscosity_dvar[_i][ph][pvar] / Utility::pow<2>(_fluid_viscosity[_i][ph]);
247 128452 : if (_multiply_by_density)
248 : {
249 128452 : const Real densityprime = (*_dfluid_density_node_dvar)[_i][ph][pvar];
250 128452 : mobprime = densityprime * mob + (*_fluid_density_node)[_i][ph] * mobprime;
251 128452 : mob *= (*_fluid_density_node)[_i][ph];
252 : }
253 128452 : if (_include_relperm)
254 : {
255 125920 : const Real relpermprime = (*_drelative_permeability_dvar)[_i][ph][pvar];
256 125920 : mobprime = relpermprime * mob + (*_relative_permeability)[_i][ph] * mobprime;
257 125920 : mob *= (*_relative_permeability)[_i][ph];
258 : }
259 :
260 128452 : const Real pre = prefactor(ph);
261 : const Real preprime =
262 128452 : (_flux_type == FluxTypeChoiceEnum::FLUID ? (*_dmass_fractions_dvar)[_i][ph][_sp][pvar]
263 37600 : : (*_denthalpy_dvar)[_i][ph][pvar]);
264 :
265 128452 : advective_term_prime -= preprime * mob * dar + pre * mobprime * dar + pre * mob * darcyprime;
266 : }
267 : }
268 766120 : if (_flux_type == FluxTypeChoiceEnum::FLUID)
269 523560 : return _test[_i][_qp] * advective_term_prime * _multiplier;
270 :
271 : const Real conduction_term_prime =
272 242560 : -_normals[_qp] *
273 242560 : ((*_dthermal_conductivity_dvar)[_qp][pvar] * (*_grad_t)[_qp] +
274 242560 : (*_thermal_conductivity)[_qp] * (*_dgrad_t_dvar)[_qp][pvar]) *
275 242560 : _phi[_j][_qp] -
276 : _normals[_qp] *
277 242560 : ((*_thermal_conductivity)[_qp] * (*_dgrad_t_dgradvar)[_qp][pvar] * _grad_phi[_j][_qp]);
278 242560 : return _test[_i][_qp] * (conduction_term_prime + advective_term_prime) * _multiplier;
279 : }
280 :
281 : Real
282 1153216 : PorousFlowOutflowBC::darcy(unsigned ph) const
283 : {
284 1153216 : return _normals[_qp] *
285 1153216 : (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
286 : }
287 :
288 : Real
289 1810112 : PorousFlowOutflowBC::mobility(unsigned ph) const
290 : {
291 1810112 : return (_multiply_by_density ? (*_fluid_density_node)[_i][ph] : 1.0) *
292 1810112 : (_include_relperm ? (*_relative_permeability)[_i][ph] : 1.0) / _fluid_viscosity[_i][ph];
293 : }
294 :
295 : Real
296 1938564 : PorousFlowOutflowBC::prefactor(unsigned ph) const
297 : {
298 1938564 : return (_flux_type == FluxTypeChoiceEnum::FLUID ? (*_mass_fractions)[_i][ph][_sp]
299 475904 : : (*_enthalpy)[_i][ph]);
300 : }
|
