www.mooseframework.org
PorousFlowSink.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "PorousFlowSink.h"
9 
10 // MOOSE includes
11 #include "MooseVariable.h"
12 
13 #include "libmesh/quadrature.h"
14 
15 // C++ includes
16 #include <iostream>
17 
18 template <>
19 InputParameters
21 {
22  InputParameters params = validParams<IntegratedBC>();
23  params.addRequiredParam<UserObjectName>(
24  "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
25  params.addParam<unsigned int>("fluid_phase",
26  "If supplied, then this BC will potentially be a function of fluid "
27  "pressure, and you can use mass_fraction_component, use_mobility, "
28  "use_relperm, use_enthalpy and use_energy. If not supplied, then "
29  "this BC can only be a function of temperature");
30  params.addParam<unsigned int>("mass_fraction_component",
31  "The index corresponding to a fluid "
32  "component. If supplied, the flux will "
33  "be multiplied by the nodal mass "
34  "fraction for the component");
35  params.addParam<bool>("use_mobility",
36  false,
37  "If true, then fluxes are multiplied by "
38  "(density*permeability_nn/viscosity), where the "
39  "'_nn' indicates the component normal to the "
40  "boundary. In this case bare_flux is measured in "
41  "Pa.m^-1. This can be used in conjunction with "
42  "other use_*");
43  params.addParam<bool>("use_relperm",
44  false,
45  "If true, then fluxes are multiplied by relative "
46  "permeability. This can be used in conjunction with "
47  "other use_*");
48  params.addParam<bool>("use_enthalpy",
49  false,
50  "If true, then fluxes are multiplied by enthalpy. "
51  "In this case bare_flux is measured in kg.m^-2.s^-1 "
52  "/ (J.kg). This can be used in conjunction with "
53  "other use_*");
54  params.addParam<bool>("use_internal_energy",
55  false,
56  "If true, then fluxes are multiplied by fluid internal energy. "
57  " In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg). "
58  " This can be used in conjunction with other use_*");
59  params.addParam<bool>("use_thermal_conductivity",
60  false,
61  "If true, then fluxes are multiplied by "
62  "thermal conductivity projected onto "
63  "the normal direction. This can be "
64  "used in conjunction with other use_*");
65  params.addParam<FunctionName>(
66  "flux_function",
67  1.0,
68  "The flux. The flux is OUT of the medium: hence positive values of "
69  "this function means this BC will act as a SINK, while negative values "
70  "indicate this flux will be a SOURCE. The functional form is useful "
71  "for spatially or temporally varying sinks. Without any use_*, this "
72  "function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case "
73  "with only heat and no fluids)");
74  params.addClassDescription("Applies a flux sink to a boundary.");
75  return params;
76 }
77 
78 PorousFlowSink::PorousFlowSink(const InputParameters & parameters)
79  : IntegratedBC(parameters),
80  _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
81  _involves_fluid(isParamValid("fluid_phase")),
82  _ph(_involves_fluid ? getParam<unsigned int>("fluid_phase") : 0),
83  _use_mass_fraction(isParamValid("mass_fraction_component")),
84  _has_mass_fraction(
85  hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
86  hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
87  "dPorousFlow_mass_frac_nodal_dvar")),
88  _sp(_use_mass_fraction ? getParam<unsigned int>("mass_fraction_component") : 0),
89  _use_mobility(getParam<bool>("use_mobility")),
90  _has_mobility(
91  hasMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp") &&
92  hasMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar") &&
93  hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
94  hasMaterialProperty<std::vector<std::vector<Real>>>(
95  "dPorousFlow_fluid_phase_density_nodal_dvar") &&
96  hasMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal") &&
97  hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
98  _use_relperm(getParam<bool>("use_relperm")),
99  _has_relperm(hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
100  hasMaterialProperty<std::vector<std::vector<Real>>>(
101  "dPorousFlow_relative_permeability_nodal_dvar")),
102  _use_enthalpy(getParam<bool>("use_enthalpy")),
103  _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
104  hasMaterialProperty<std::vector<std::vector<Real>>>(
105  "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
106  _use_internal_energy(getParam<bool>("use_internal_energy")),
107  _has_internal_energy(
108  hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_internal_energy_nodal") &&
109  hasMaterialProperty<std::vector<std::vector<Real>>>(
110  "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")),
111  _use_thermal_conductivity(getParam<bool>("use_thermal_conductivity")),
112  _has_thermal_conductivity(
113  hasMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp") &&
114  hasMaterialProperty<std::vector<RealTensorValue>>(
115  "dPorousFlow_thermal_conductivity_qp_dvar")),
116  _m_func(getFunction("flux_function")),
117  _permeability(_has_mobility
118  ? &getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")
119  : nullptr),
120  _dpermeability_dvar(
121  _has_mobility
122  ? &getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")
123  : nullptr),
124  _dpermeability_dgradvar(_has_mobility
125  ? &getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
126  "dPorousFlow_permeability_qp_dgradvar")
127  : nullptr),
128  _fluid_density_node(
129  _has_mobility
130  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")
131  : nullptr),
132  _dfluid_density_node_dvar(_has_mobility
133  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
134  "dPorousFlow_fluid_phase_density_nodal_dvar")
135  : nullptr),
136  _fluid_viscosity(_has_mobility
137  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
138  : nullptr),
139  _dfluid_viscosity_dvar(_has_mobility
140  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
141  "dPorousFlow_viscosity_nodal_dvar")
142  : nullptr),
143  _relative_permeability(
144  _has_relperm
145  ? &getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal")
146  : nullptr),
147  _drelative_permeability_dvar(_has_relperm
148  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
149  "dPorousFlow_relative_permeability_nodal_dvar")
150  : nullptr),
151  _mass_fractions(
152  _has_mass_fraction
153  ? &getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")
154  : nullptr),
155  _dmass_fractions_dvar(_has_mass_fraction
156  ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
157  "dPorousFlow_mass_frac_nodal_dvar")
158  : nullptr),
159  _enthalpy(
160  _has_enthalpy
161  ? &getMaterialPropertyByName<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal")
162  : nullptr),
163  _denthalpy_dvar(_has_enthalpy
164  ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
165  "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
166  : nullptr),
167  _internal_energy(_has_internal_energy
168  ? &getMaterialPropertyByName<std::vector<Real>>(
169  "PorousFlow_fluid_phase_internal_energy_nodal")
170  : nullptr),
171  _dinternal_energy_dvar(_has_internal_energy
172  ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
173  "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
174  : nullptr),
175  _thermal_conductivity(
176  _has_thermal_conductivity
177  ? &getMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp")
178  : nullptr),
179  _dthermal_conductivity_dvar(_has_thermal_conductivity
180  ? &getMaterialProperty<std::vector<RealTensorValue>>(
181  "dPorousFlow_thermal_conductivity_qp_dvar")
182  : nullptr)
183 {
185  mooseError("PorousFlowSink: The Dictator declares that the number of fluid phases is ",
187  ", but you have set the fluid_phase to ",
188  _ph,
189  ". You must try harder.");
192  mooseError("PorousFlowSink: To use_mass_fraction, use_mobility, use_relperm, use_enthalpy or "
193  "use_internal_energy, you must provide a fluid phase number");
195  mooseError("PorousFlowSink: The Dictator declares that the number of fluid components is ",
197  ", but you have set the mass_fraction_component to ",
198  _sp,
199  ". Please be assured that the Dictator has noted your error.");
201  mooseError("PorousFlowSink: You have used the use_mass_fraction flag, but you have no "
202  "mass_fraction Material");
204  mooseError("PorousFlowSink: You have used the use_mobility flag, but there are not the "
205  "required Materials for this");
206  if (_use_relperm && !_has_relperm)
207  mooseError(
208  "PorousFlowSink: You have used the use_relperm flag, but you have no relperm Material");
210  mooseError(
211  "PorousFlowSink: You have used the use_enthalpy flag, but you have no enthalpy Material");
213  mooseError("PorousFlowSink: You have used the use_internal_energy flag, but you have no "
214  "internal_energy Material");
216  mooseError("PorousFlowSink: You have used the use_thermal_conductivity flag, but you have no "
217  "thermal_conductivity Material");
218 }
219 
220 Real
222 {
223  Real flux = _test[_i][_qp] * multiplier();
224  if (_use_mobility)
225  {
226  const Real k =
227  ((*_permeability)[_qp] * _normals[_qp]) * _normals[_qp]; // do not upwind permeability
228  flux *= (*_fluid_density_node)[_i][_ph] * k / (*_fluid_viscosity)[_i][_ph];
229  }
230  if (_use_relperm)
231  flux *= (*_relative_permeability)[_i][_ph];
232  if (_use_mass_fraction)
233  flux *= (*_mass_fractions)[_i][_ph][_sp];
234  if (_use_enthalpy)
235  flux *= (*_enthalpy)[_i][_ph];
237  flux *= (*_internal_energy)[_i][_ph];
239  flux *= ((*_thermal_conductivity)[_qp] * _normals[_qp]) *
240  _normals[_qp]; // do not upwind thermal_conductivity
241 
242  return flux;
243 }
244 
245 Real
247 {
248  return jac(_var.number());
249 }
250 
251 Real
253 {
254  return jac(jvar);
255 }
256 
257 Real
258 PorousFlowSink::jac(unsigned int jvar) const
259 {
261  return 0.0;
262  const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
263 
264  // For _i != _j, note:
265  // since the only non-upwinded contribution to the residual is
266  // from the permeability and thermal_conductivity, the only contribution
267  // of the residual at node _i from changing jvar at node _j is through
268  // the derivative of permeability or thermal_conductivity
269 
270  Real flux = _test[_i][_qp] * multiplier();
271  Real deriv = _test[_i][_qp] * (_i != _j ? 0.0 : dmultiplier_dvar(pvar));
272 
273  if (_use_mobility)
274  {
275  const Real k = ((*_permeability)[_qp] * _normals[_qp]) * _normals[_qp];
276  const Real mob = (*_fluid_density_node)[_i][_ph] * k / (*_fluid_viscosity)[_i][_ph];
277  RealTensorValue ktprime = (*_dpermeability_dvar)[_qp][pvar] * _phi[_j][_qp];
278  for (unsigned i = 0; i < LIBMESH_DIM; ++i)
279  ktprime += (*_dpermeability_dgradvar)[_qp][i][pvar] * _grad_phi[_j][_qp](i);
280  const Real kprime = (ktprime * _normals[_qp]) * _normals[_qp];
281 
282  Real mobprime = (*_fluid_density_node)[_i][_ph] * kprime / (*_fluid_viscosity)[_i][_ph];
283  mobprime +=
284  (_i != _j
285  ? 0.0
286  : (*_dfluid_density_node_dvar)[_i][_ph][pvar] * k / (*_fluid_viscosity)[_i][_ph] -
287  (*_fluid_density_node)[_i][_ph] * k * (*_dfluid_viscosity_dvar)[_i][_ph][pvar] /
288  std::pow((*_fluid_viscosity)[_i][_ph], 2));
289  deriv = mob * deriv + mobprime * flux;
290  flux *= mob;
291  }
292  if (_use_relperm)
293  {
294  const Real relperm_prime = (_i != _j ? 0.0 : (*_drelative_permeability_dvar)[_i][_ph][pvar]);
295  deriv = (*_relative_permeability)[_i][_ph] * deriv + relperm_prime * flux;
296  flux *= (*_relative_permeability)[_i][_ph];
297  }
298  if (_use_mass_fraction)
299  {
300  const Real mf_prime = (_i != _j ? 0.0 : (*_dmass_fractions_dvar)[_i][_ph][_sp][pvar]);
301  deriv = (*_mass_fractions)[_i][_ph][_sp] * deriv + mf_prime * flux;
302  flux *= (*_mass_fractions)[_i][_ph][_sp];
303  }
304  if (_use_enthalpy)
305  {
306  const Real en_prime = (_i != _j ? 0.0 : (*_denthalpy_dvar)[_i][_ph][pvar]);
307  deriv = (*_enthalpy)[_i][_ph] * deriv + en_prime * flux;
308  flux *= (*_enthalpy)[_i][_ph];
309  }
311  {
312  const Real ie_prime = (_i != _j ? 0.0 : (*_dinternal_energy_dvar)[_i][_ph][pvar]);
313  deriv = (*_internal_energy)[_i][_ph] * deriv + ie_prime * flux;
314  flux *= (*_internal_energy)[_i][_ph];
315  }
317  {
318  const Real tc = ((*_thermal_conductivity)[_qp] * _normals[_qp]) * _normals[_qp];
319  const RealTensorValue tctprime = (*_dthermal_conductivity_dvar)[_qp][pvar] * _phi[_j][_qp];
320  const Real tcprime = (tctprime * _normals[_qp]) * _normals[_qp];
321  deriv = tc * deriv + tcprime * flux;
322  // don't need this: flux *= tc;
323  }
324  return deriv;
325 }
326 
327 Real
329 {
330  return _m_func.value(_t, _q_point[_qp]);
331 }
332 
333 Real
334 PorousFlowSink::dmultiplier_dvar(unsigned int /*pvar*/) const
335 {
336  return 0.0;
337 }
const bool _has_enthalpy
whether there is an "enthalpy" Material. This is just for error checking
const bool _involves_fluid
Whether this BC involves fluid (whether the user has supplied a fluid phase number) ...
const unsigned int _sp
The component number (only used if _use_mass_fraction==true)
const unsigned int _ph
The phase number.
const MaterialProperty< std::vector< std::vector< RealTensorValue > > > *const _dpermeability_dgradvar
d(Permeability)/d(grad(PorousFlow variable))
const PorousFlowDictator & _dictator
PorousFlow UserObject.
const bool _has_internal_energy
whether there is an "internal_energy" Material. This is just for error checking
const MaterialProperty< std::vector< Real > > *const _fluid_viscosity
Viscosity of each component in each phase.
const bool _has_thermal_conductivity
whether there is an "thermal_conductivity" Material. This is just for error checking ...
const bool _has_mobility
Whether there are Materials that can form "mobility". This is just for error checking.
const bool _use_mass_fraction
Whether the flux will be multiplied by the mass fraction.
unsigned int numComponents() const
the number of fluid components
InputParameters validParams< PorousFlowSink >()
const bool _has_relperm
Whether there is a "relperm" Material. This is just for error checking.
const bool _use_relperm
whether to multiply the sink flux by relative permeability
const bool _use_mobility
whether to multiply the sink flux by permeability*density/viscosity
virtual Real computeQpJacobian() override
Real jac(unsigned int jvar) const
derivative of residual with respect to the jvar variable
const bool _use_enthalpy
whether to multiply the sink flux by enthalpy
Function & _m_func
The flux.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
virtual Real dmultiplier_dvar(unsigned int pvar) const
d(multiplier)/d(Porous flow variable pvar)
This holds maps between the nonlinear variables used in a PorousFlow simulation and the variable numb...
unsigned int numPhases() const
the number of fluid phases
const bool _has_mass_fraction
Whether there is a "mass_fraction" Material. This is just for error checking.
const bool _use_internal_energy
whether to multiply the sink flux by internal_energy
virtual Real computeQpResidual() override
const bool _use_thermal_conductivity
whether to multiply the sink flux by thermal_conductivity
bool notPorousFlowVariable(unsigned int moose_var_num) const
returns true if moose_var_num is not a porous flow variabe
PorousFlowSink(const InputParameters &parameters)
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
the PorousFlow variable number
virtual Real multiplier() const
The flux gets multiplied by this quantity.