www.mooseframework.org
Q2PBorehole.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 "Q2PBorehole.h"
9 #include "RotationMatrix.h"
10 
11 template <>
12 InputParameters
14 {
15  InputParameters params = validParams<PeacemanBorehole>();
16  params.addRequiredParam<UserObjectName>(
17  "fluid_density",
18  "A RichardsDensity UserObject that defines the fluid density as a function of pressure.");
19  params.addRequiredParam<UserObjectName>(
20  "fluid_relperm",
21  "A RichardsRelPerm UserObject (eg RichardsRelPermPower) that defines the "
22  "fluid relative permeability as a function of the saturation Variable.");
23  params.addRequiredCoupledVar("other_var",
24  "The other variable in the 2-phase system. If "
25  "Variable=porepressure, the other_var=saturation, and "
26  "vice-versa.");
27  params.addRequiredParam<bool>("var_is_porepressure",
28  "This flag is needed to correctly calculate the Jacobian entries. "
29  "If set to true, this Sink will extract fluid from the phase with "
30  "porepressure as its Variable (usually the liquid phase). If set "
31  "to false, this Sink will extract fluid from the phase with "
32  "saturation as its variable (usually the gas phase)");
33  params.addRequiredParam<Real>("fluid_viscosity", "The fluid dynamic viscosity");
34  params.addClassDescription("Approximates a borehole in the mesh with given bottomhole pressure, "
35  "and radii using a number of point sinks whose positions are read "
36  "from a file. This DiracKernel is for use by Q2P models");
37  return params;
38 }
39 
40 Q2PBorehole::Q2PBorehole(const InputParameters & parameters)
41  : PeacemanBorehole(parameters),
42  _density(getUserObject<RichardsDensity>("fluid_density")),
43  _relperm(getUserObject<RichardsRelPerm>("fluid_relperm")),
44  _other_var_nodal(coupledNodalValue("other_var")),
45  _other_var_num(coupled("other_var")),
46  _var_is_pp(getParam<bool>("var_is_porepressure")),
47  _viscosity(getParam<Real>("fluid_viscosity")),
48  _permeability(getMaterialProperty<RealTensorValue>("permeability")),
49  _num_nodes(0),
50  _pp(0),
51  _sat(0),
52  _mobility(0),
53  _dmobility_dp(0),
54  _dmobility_ds(0)
55 {
56 }
57 
58 void
60 {
62 
63  // set _pp and _sat variables
64  _pp.resize(_num_nodes);
65  _sat.resize(_num_nodes);
66  if (_var_is_pp)
67  {
68  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
69  {
70  _pp[nodenum] = _var.nodalSln()[nodenum];
71  _sat[nodenum] = _other_var_nodal[nodenum];
72  }
73  }
74  else
75  {
76  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
77  {
78  _pp[nodenum] = _other_var_nodal[nodenum];
79  _sat[nodenum] = _var.nodalSln()[nodenum];
80  }
81  }
82 
83  Real density;
84  Real ddensity_dp;
85  Real relperm;
86  Real drelperm_ds;
87  _mobility.resize(_num_nodes);
88  _dmobility_dp.resize(_num_nodes);
89  _dmobility_ds.resize(_num_nodes);
90  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
91  {
92  density = _density.density(_pp[nodenum]);
93  ddensity_dp = _density.ddensity(_pp[nodenum]);
94  relperm = _relperm.relperm(_sat[nodenum]);
95  drelperm_ds = _relperm.drelperm(_sat[nodenum]);
96  _mobility[nodenum] = density * relperm / _viscosity;
97  _dmobility_dp[nodenum] = ddensity_dp * relperm / _viscosity;
98  _dmobility_ds[nodenum] = density * drelperm_ds / _viscosity;
99  }
100 }
101 
102 void
104 {
106  DiracKernel::computeResidual();
107 }
108 
109 Real
111 {
112  const Real character = _character.value(_t, _q_point[_qp]);
113  if (character == 0.0)
114  return 0.0;
115 
116  const Real bh_pressure =
117  _p_bot +
118  _unit_weight *
119  (_q_point[_qp] -
120  _bottom_point); // really want to use _q_point instaed of _current_point, i think?!
121 
122  // Get the ID we initially assigned to this point
123  const unsigned current_dirac_ptid = currentPointCachedID();
124 
125  // If getting the ID failed, fall back to the old bodge!
126  // if (current_dirac_ptid == libMesh::invalid_uint)
127  // current_dirac_ptid = (_zs.size() > 2) ? 1 : 0;
128 
129  Real outflow(0.0); // this is the flow rate from porespace out of the system
130 
131  Real wc(0.0);
132  if (current_dirac_ptid > 0)
133  // contribution from half-segment "behind" this point (must have >1 point for
134  // current_dirac_ptid>0)
135  {
136  wc = wellConstant(_permeability[0],
137  _rot_matrix[current_dirac_ptid - 1],
138  _half_seg_len[current_dirac_ptid - 1],
139  _current_elem,
140  _rs[current_dirac_ptid]);
141  if ((character < 0.0 && _pp[_i] < bh_pressure) || (character > 0.0 && _pp[_i] > bh_pressure))
142  // injection, so outflow<0 || // production, so outflow>0
143  outflow +=
144  _test[_i][_qp] * std::abs(character) * wc * _mobility[_i] * (_pp[_i] - bh_pressure);
145  }
146 
147  if (current_dirac_ptid + 1 < _zs.size() || _zs.size() == 1)
148  // contribution from half-segment "ahead of" this point, or we only have one point
149  {
150  wc = wellConstant(_permeability[0],
151  _rot_matrix[current_dirac_ptid],
152  _half_seg_len[current_dirac_ptid],
153  _current_elem,
154  _rs[current_dirac_ptid]);
155  if ((character < 0.0 && _pp[_i] < bh_pressure) || (character > 0.0 && _pp[_i] > bh_pressure))
156  // injection, so outflow<0 || // production, so outflow>0
157  outflow +=
158  _test[_i][_qp] * std::abs(character) * wc * _mobility[_i] * (_pp[_i] - bh_pressure);
159  }
160 
162  outflow * _dt); // this is not thread safe, but DiracKernel's aren't currently threaded
163  return outflow;
164 }
165 
166 void
168 {
170  DiracKernel::computeJacobian();
171 }
172 
173 Real
175 {
176  return jac(_var.number());
177 }
178 
179 Real
181 {
182  if (jvar == _other_var_num || jvar == _var.number())
183  return jac(jvar);
184  return 0.0;
185 }
186 
187 Real
188 Q2PBorehole::jac(unsigned int jvar)
189 {
190  if (_i != _j)
191  return 0.0;
192 
193  const Real character = _character.value(_t, _q_point[_qp]);
194  if (character == 0.0)
195  return 0.0;
196 
197  const Real bh_pressure =
198  _p_bot +
199  _unit_weight *
200  (_q_point[_qp] -
201  _bottom_point); // really want to use _q_point instaed of _current_point, i think?!
202 
203  const Real phi = 1;
204 
205  // Get the ID we initially assigned to this point
206  const unsigned current_dirac_ptid = currentPointCachedID();
207 
208  // If getting the ID failed, fall back to the old bodge!
209  // if (current_dirac_ptid == libMesh::invalid_uint)
210  // current_dirac_ptid = (_zs.size() > 2) ? 1 : 0;
211 
212  Real outflowp(0.0);
213 
214  const bool deriv_wrt_pp =
215  (_var_is_pp && (jvar == _var.number())) || (!_var_is_pp && (jvar == _other_var_num));
216 
217  Real wc(0.0);
218  if (current_dirac_ptid > 0)
219  // contribution from half-segment "behind" this point
220  {
221  wc = wellConstant(_permeability[0],
222  _rot_matrix[current_dirac_ptid - 1],
223  _half_seg_len[current_dirac_ptid - 1],
224  _current_elem,
225  _rs[current_dirac_ptid]);
226  if ((character < 0.0 && _pp[_i] < bh_pressure) || (character > 0.0 && _pp[_i] > bh_pressure))
227  {
228  // injection, so outflow<0 || // production, so outflow>0
229  if (deriv_wrt_pp)
230  outflowp += std::abs(character) * wc *
231  (_mobility[_i] * phi + _dmobility_dp[_i] * phi * (_pp[_i] - bh_pressure));
232  else
233  outflowp += std::abs(character) * wc * _dmobility_ds[_i] * phi * (_pp[_i] - bh_pressure);
234  }
235  }
236 
237  if (current_dirac_ptid < _zs.size() - 1 || _zs.size() == 1)
238  // contribution from half-segment "ahead of" this point
239  {
240  wc = wellConstant(_permeability[0],
241  _rot_matrix[current_dirac_ptid],
242  _half_seg_len[current_dirac_ptid],
243  _current_elem,
244  _rs[current_dirac_ptid]);
245  if ((character < 0.0 && _pp[_i] < bh_pressure) || (character > 0.0 && _pp[_i] > bh_pressure))
246  {
247  // injection, so outflow<0 || // production, so outflow>0
248  if (deriv_wrt_pp)
249  outflowp += std::abs(character) * wc *
250  (_mobility[_i] * phi + _dmobility_dp[_i] * phi * (_pp[_i] - bh_pressure));
251  else
252  outflowp += std::abs(character) * wc * _dmobility_ds[_i] * phi * (_pp[_i] - bh_pressure);
253  }
254  }
255 
256  return _test[_i][_qp] * outflowp;
257 }
unsigned int _num_nodes
number of nodes in this element.
Definition: Q2PBorehole.h:91
InputParameters validParams< PeacemanBorehole >()
const RichardsRelPerm & _relperm
fluid relative permeability
Definition: Q2PBorehole.h:73
virtual Real drelperm(Real seff) const =0
derivative of relative permeability wrt effective saturation This must be over-ridden in your derived...
std::vector< Real > _dmobility_dp
nodal d(mobility)/d(porepressure)
Definition: Q2PBorehole.h:103
Real wellConstant(const RealTensorValue &perm, const RealTensorValue &rot, const Real &half_len, const Elem *ele, const Real &rad)
Calculates Peaceman&#39;s form of the borehole well constant Z Chen, Y Zhang, Well flow models for variou...
Base class for Richards relative permeability classes that provide relative permeability as a functio...
const unsigned int _other_var_num
the variable number of the other variable
Definition: Q2PBorehole.h:79
std::vector< Real > _rs
radii of the borehole
virtual Real computeQpResidual()
Computes the Qp residual.
Definition: Q2PBorehole.C:110
const std::string density
Definition: NS.h:15
const VariableValue & _other_var_nodal
the other variable in the 2-phase system (this is saturation if Variable=porepressure, and viceversa)
Definition: Q2PBorehole.h:76
Approximates a borehole by a sequence of Dirac Points.
std::vector< Real > _mobility
nodal mobility
Definition: Q2PBorehole.h:100
Q2PBorehole(const InputParameters &parameters)
Creates a new Q2PBorehole This sets all the variables, but also reads the file containing the lines o...
Definition: Q2PBorehole.C:40
virtual void computeResidual()
Computes the residual.
Definition: Q2PBorehole.C:103
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...
void prepareNodalValues()
calculates the nodal values of pressure, mobility, and derivatives thereof
Definition: Q2PBorehole.C:59
InputParameters validParams< Q2PBorehole >()
Definition: Q2PBorehole.C:13
virtual void computeJacobian()
Computes the Jacobian.
Definition: Q2PBorehole.C:167
const Real _viscosity
viscosity
Definition: Q2PBorehole.h:85
std::vector< Real > _sat
nodal saturation
Definition: Q2PBorehole.h:97
const RealVectorValue _unit_weight
unit weight of fluid in borehole (for calculating bottomhole pressure at each Dirac Point) ...
virtual Real computeQpJacobian()
Computes the diagonal part of the jacobian.
Definition: Q2PBorehole.C:174
void add(Real contrib)
adds contrib to _total
RichardsSumQuantity & _total_outflow_mass
This is used to hold the total fluid flowing into the borehole Hence, it is positive for production w...
const RichardsDensity & _density
fluid density
Definition: Q2PBorehole.h:70
virtual Real ddensity(Real p) const =0
derivative of fluid density wrt porepressure This must be over-ridden in derived classes to provide a...
const Real _p_bot
bottomhole pressure of borehole
std::vector< Real > _dmobility_ds
nodal d(mobility)/d(saturation)
Definition: Q2PBorehole.h:106
const bool _var_is_pp
whether the Variable for this BC is porepressure or not
Definition: Q2PBorehole.h:82
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
Definition: Q2PBorehole.h:88
Real jac(unsigned int jvar)
Calculates Jacobian.
Definition: Q2PBorehole.C:188
Point _bottom_point
the bottom point of the borehole (where bottom_pressure is defined)
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
std::vector< RealTensorValue > _rot_matrix
rotation matrix used in well_constant calculation
std::vector< Real > _half_seg_len
0.5*(length of polyline segments between points)
std::vector< Real > _zs
z points of borehole
std::vector< Real > _pp
nodal porepressure
Definition: Q2PBorehole.h:94
virtual Real computeQpOffDiagJacobian(unsigned int jvar)
Computes the off-diagonal part of the jacobian.
Definition: Q2PBorehole.C:180
Function & _character
If positive then the borehole acts as a sink (producion well) for porepressure > borehole pressure...