www.mooseframework.org
Q2PPiecewiseLinearSink.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 
9 
10 // MOOSE includes
11 #include "MooseVariable.h"
12 
13 // C++ includes
14 #include <iostream>
15 
16 template <>
17 InputParameters
19 {
20  InputParameters params = validParams<IntegratedBC>();
21  params.addRequiredParam<bool>(
22  "use_mobility",
23  "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), "
24  "where the '_nn' indicates the component normal to the boundary. In this "
25  "case bare_flux is measured in Pa.s^-1. This can be used in conjunction "
26  "with use_relperm.");
27  params.addRequiredParam<bool>("use_relperm",
28  "If true, then fluxes are multiplied by relative "
29  "permeability. This can be used in conjunction "
30  "with use_mobility");
31  params.addRequiredParam<std::vector<Real>>(
32  "pressures", "Tuple of pressure values. Must be monotonically increasing.");
33  params.addRequiredParam<std::vector<Real>>(
34  "bare_fluxes",
35  "Tuple of flux values (measured in kg.m^-2.s^-1 for use_mobility=false, and "
36  "in Pa.s^-1 if use_mobility=true). This flux is OUT of the medium: hence "
37  "positive values of flux means this will be a SINK, while negative values "
38  "indicate this flux will be a SOURCE. A piecewise-linear fit is performed to "
39  "the (pressure,bare_fluxes) pairs to obtain the flux at any arbitrary "
40  "pressure, and the first or last bare_flux values are used if the quad-point "
41  "pressure falls outside this range.");
42  params.addParam<FunctionName>("multiplying_fcn",
43  1.0,
44  "If this function is provided, the flux "
45  "will be multiplied by this function. "
46  "This is useful for spatially or "
47  "temporally varying sinks");
48  params.addRequiredParam<UserObjectName>(
49  "fluid_density",
50  "A RichardsDensity UserObject that defines the fluid density as a function of pressure.");
51  params.addRequiredParam<UserObjectName>(
52  "fluid_relperm",
53  "A RichardsRelPerm UserObject (eg RichardsRelPermPower) that defines the "
54  "fluid relative permeability as a function of the saturation Variable.");
55  params.addRequiredCoupledVar("other_var",
56  "The other variable in the 2-phase system. If "
57  "Variable=porepressure, the other_var=saturation, and "
58  "vice-versa.");
59  params.addRequiredParam<bool>("var_is_porepressure",
60  "This flag is needed to correctly calculate the Jacobian entries. "
61  "If set to true, this Sink will extract fluid from the phase with "
62  "porepressure as its Variable (usually the liquid phase). If set "
63  "to false, this Sink will extract fluid from the phase with "
64  "saturation as its variable (usually the gas phase)");
65  params.addRequiredParam<Real>("fluid_viscosity", "The fluid dynamic viscosity");
66  params.addClassDescription("Sink of fluid, controlled by (pressure, bare_fluxes) interpolation. "
67  "This is for use in Q2P models");
68  return params;
69 }
70 
71 Q2PPiecewiseLinearSink::Q2PPiecewiseLinearSink(const InputParameters & parameters)
72  : IntegratedBC(parameters),
73  _use_mobility(getParam<bool>("use_mobility")),
74  _use_relperm(getParam<bool>("use_relperm")),
75  _sink_func(getParam<std::vector<Real>>("pressures"),
76  getParam<std::vector<Real>>("bare_fluxes")),
77  _m_func(getFunction("multiplying_fcn")),
78  _density(getUserObject<RichardsDensity>("fluid_density")),
79  _relperm(getUserObject<RichardsRelPerm>("fluid_relperm")),
80  _other_var_nodal(coupledNodalValue("other_var")),
81  _other_var_num(coupled("other_var")),
82  _var_is_pp(getParam<bool>("var_is_porepressure")),
83  _viscosity(getParam<Real>("fluid_viscosity")),
84  _permeability(getMaterialProperty<RealTensorValue>("permeability")),
85  _num_nodes(0),
86  _pp(0),
87  _sat(0),
88  _nodal_density(0),
89  _dnodal_density_dp(0),
90  _nodal_relperm(0),
91  _dnodal_relperm_ds(0)
92 {
93 }
94 
95 void
97 {
99 
100  // set _pp and _sat variables
101  _pp.resize(_num_nodes);
102  _sat.resize(_num_nodes);
103  if (_var_is_pp)
104  {
105  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
106  {
107  _pp[nodenum] = _var.nodalSln()[nodenum];
108  _sat[nodenum] = _other_var_nodal[nodenum];
109  }
110  }
111  else
112  {
113  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
114  {
115  _pp[nodenum] = _other_var_nodal[nodenum];
116  _sat[nodenum] = _var.nodalSln()[nodenum];
117  }
118  }
119 
120  // calculate derived things
121  if (_use_mobility)
122  {
123  _nodal_density.resize(_num_nodes);
125  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
126  {
127  _nodal_density[nodenum] = _density.density(_pp[nodenum]);
128  _dnodal_density_dp[nodenum] = _density.ddensity(_pp[nodenum]);
129  }
130  }
131 
132  if (_use_relperm)
133  {
134  _nodal_relperm.resize(_num_nodes);
136  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
137  {
138  _nodal_relperm[nodenum] = _relperm.relperm(_sat[nodenum]);
139  _dnodal_relperm_ds[nodenum] = _relperm.drelperm(_sat[nodenum]);
140  }
141  }
142 }
143 
144 void
146 {
148  IntegratedBC::computeResidual();
149 }
150 
151 Real
153 {
154  Real flux = 0;
155  Real k = 0;
156 
157  flux = _test[_i][_qp] * _sink_func.sample(_pp[_i]);
158  if (_use_mobility)
159  {
160  k = (_permeability[0] * _normals[_qp]) * _normals[_qp]; // assume that _permeability is constant
161  // throughout element so doesn't need to
162  // be upwinded
163  flux *= _nodal_density[_i] * k / _viscosity;
164  }
165  if (_use_relperm)
166  flux *= _nodal_relperm[_i];
167 
168  flux *= _m_func.value(_t, _q_point[_qp]);
169 
170  return flux;
171 }
172 
173 void
175 {
177  IntegratedBC::computeJacobian();
178 }
179 
180 Real
182 {
183  return jac(_var.number());
184 }
185 
186 void
188 {
190  IntegratedBC::computeJacobianBlock(jvar);
191 }
192 
193 Real
195 {
196  if (jvar == _var.number() || jvar == _other_var_num)
197  return jac(jvar);
198 
199  return 0.0;
200 }
201 
202 Real
203 Q2PPiecewiseLinearSink::jac(unsigned int wrt_num)
204 {
205  Real flux = 0;
206  Real deriv = 0;
207  Real k = 0;
208  Real mob = 0;
209  Real mobp = 0;
210  Real phi = 1;
211 
212  if (_i != _j)
213  return 0.0; // residual at node _i only depends on variables at that node
214 
215  flux = _sink_func.sample(_pp[_i]);
216 
217  if (_var_is_pp)
218  {
219  // derivative of the _sink_func
220  if (wrt_num == _var.number())
221  deriv = _sink_func.sampleDerivative(_pp[_i]);
222  else
223  deriv = 0;
224 
225  // add derivative of the mobility
226  if (_use_mobility)
227  {
228  k = (_permeability[0] * _normals[_qp]) * _normals[_qp];
229  mob = _nodal_density[_i] * k / _viscosity;
230  if (wrt_num == _var.number())
231  mobp = _dnodal_density_dp[_i] * k / _viscosity; // else mobp = 0
232  deriv = mob * deriv + mobp * flux;
233  flux *= mob;
234  }
235 
236  // add derivative of the relperm
237  if (_use_relperm)
238  {
239  if (wrt_num == _other_var_num)
240  deriv = _nodal_relperm[_i] * deriv + _dnodal_relperm_ds[_i] * flux;
241  else
242  deriv = _nodal_relperm[_i] * deriv;
243  }
244  }
245  else
246  {
247  // derivative of the _sink_func
248  if (wrt_num == _other_var_num)
249  deriv = _sink_func.sampleDerivative(_pp[_i]);
250  else
251  deriv = 0;
252 
253  // add derivative of the mobility
254  if (_use_mobility)
255  {
256  k = (_permeability[0] * _normals[_qp]) * _normals[_qp];
257  mob = _nodal_density[_i] * k / _viscosity;
258  if (wrt_num == _other_var_num)
259  mobp = _dnodal_density_dp[_i] * k / _viscosity; // else mobp = 0
260  deriv = mob * deriv + mobp * flux;
261  flux *= mob;
262  }
263 
264  // add derivative of the relperm
265  if (_use_relperm)
266  {
267  if (wrt_num == _var.number())
268  deriv = _nodal_relperm[_i] * deriv + _dnodal_relperm_ds[_i] * flux;
269  else
270  deriv = _nodal_relperm[_i] * deriv;
271  }
272  }
273 
274  deriv *= _m_func.value(_t, _q_point[_qp]);
275 
276  return _test[_i][_qp] * deriv * phi;
277 }
void prepareNodalValues()
calculates the nodal values of pressure, mobility, and derivatives thereof
Function & _m_func
sink flux gets multiplied by this function
std::vector< Real > _nodal_density
nodal values of fluid density
virtual Real drelperm(Real seff) const =0
derivative of relative permeability wrt effective saturation This must be over-ridden in your derived...
LinearInterpolation _sink_func
piecewise-linear function of porepressure (this defines the strength of the sink) ...
bool _use_mobility
whether to multiply the sink flux by permeability*density/viscosity
Q2PPiecewiseLinearSink(const InputParameters &parameters)
Base class for Richards relative permeability classes that provide relative permeability as a functio...
virtual void computeJacobianBlock(unsigned int jvar)
virtual Real computeQpOffDiagJacobian(unsigned int jvar)
Real jac(unsigned int wrt_num)
derivative of residual wrt the wrt_num variable
unsigned int _other_var_num
the variable number of the other variable
virtual Real density(Real p) const =0
fluid density as a function of porepressure This must be over-ridden in derived classes to provide an...
bool _var_is_pp
whether the Variable for this BC is porepressure or not
std::vector< Real > _dnodal_density_dp
d(_nodal_density)/d(porepressure)
std::vector< Real > _pp
nodal values of porepressure
const RichardsRelPerm & _relperm
fluid relative permeability
virtual Real ddensity(Real p) const =0
derivative of fluid density wrt porepressure This must be over-ridden in derived classes to provide a...
virtual Real relperm(Real seff) const =0
relative permeability as a function of effective saturation This must be over-ridden in your derived ...
const MaterialProperty< RealTensorValue > & _permeability
permeability
bool _use_relperm
whether to multiply the sink flux by relative permeability
std::vector< Real > _dnodal_relperm_ds
d(_nodal_relperm)/d(saturation)
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
InputParameters validParams< Q2PPiecewiseLinearSink >()
unsigned int _num_nodes
number of nodes in this element.
std::vector< Real > _nodal_relperm
nodal values of relative permeability
std::vector< Real > _sat
nodal values of saturation
const RichardsDensity & _density
fluid density
const VariableValue & _other_var_nodal
the other variable in the 2-phase system (this is saturation if Variable=porepressure, and viceversa)