www.mooseframework.org
PorousFlowPeacemanBorehole.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
11 #include "RotationMatrix.h"
12 #include "Function.h"
13 
15 
18 {
20  params.addRequiredParam<FunctionName>(
21  "character",
22  "If zero then borehole does nothing. If positive the borehole acts as a sink "
23  "(production well) for porepressure > borehole pressure, and does nothing "
24  "otherwise. If negative the borehole acts as a source (injection well) for "
25  "porepressure < borehole pressure, and does nothing otherwise. The flow rate "
26  "to/from the borehole is multiplied by |character|, so usually character = +/- "
27  "1, but you can specify other quantities to provide an overall scaling to the "
28  "flow if you like.");
29  params.addRequiredParam<FunctionName>("bottom_p_or_t",
30  "For function_of=pressure, this function is the "
31  "pressure at the bottom of the borehole, "
32  "otherwise it is the temperature at the bottom of "
33  "the borehole.");
35  "unit_weight",
36  "(fluid_density*gravitational_acceleration) as a vector pointing downwards. "
37  "Note that the borehole pressure at a given z position is bottom_p_or_t + "
38  "unit_weight*(q - q_bottom), where q=(x,y,z) and q_bottom=(x,y,z) of the "
39  "bottom point of the borehole. The analogous formula holds for "
40  "function_of=temperature. If you don't want bottomhole pressure (or "
41  "temperature) to vary in the borehole just set unit_weight=0. Typical value "
42  "is = (0,0,-1E4), for water");
43  params.addParam<Real>("re_constant",
44  0.28,
45  "The dimensionless constant used in evaluating the borehole effective "
46  "radius. This depends on the meshing scheme. Peacemann "
47  "finite-difference calculations give 0.28, while for rectangular finite "
48  "elements the result is closer to 0.1594. (See Eqn(4.13) of Z Chen, Y "
49  "Zhang, Well flow models for various numerical methods, Int J Num "
50  "Analysis and Modeling, 3 (2008) 375-388.)");
51  params.addParam<Real>("well_constant",
52  -1.0,
53  "Usually this is calculated internally from the element geometry, the "
54  "local borehole direction and segment length, and the permeability. "
55  "However, if this parameter is given as a positive number then this "
56  "number is used instead of the internal calculation. This speeds up "
57  "computation marginally. re_constant becomes irrelevant");
58  params.addClassDescription(
59  "Approximates a borehole in the mesh using the Peaceman approach, ie "
60  "using a number of point sinks with given radii whose positions are "
61  "read from a file. NOTE: if you are using PorousFlowPorosity that depends on volumetric "
62  "strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal "
63  "Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be "
64  "computed");
65  return params;
66 }
67 
69  : PorousFlowLineSink(parameters),
70  _character(getFunction("character")),
71  _p_bot(getFunction("bottom_p_or_t")),
72  _unit_weight(getParam<RealVectorValue>("unit_weight")),
73  _re_constant(getParam<Real>("re_constant")),
74  _well_constant(getParam<Real>("well_constant")),
75  _has_permeability(
76  hasMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp") &&
77  hasMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
78  _has_thermal_conductivity(
79  hasMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp") &&
80  hasMaterialProperty<std::vector<RealTensorValue>>(
81  "dPorousFlow_thermal_conductivity_qp_dvar")),
82  _perm_or_cond(_p_or_t == PorTchoice::pressure
83  ? getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")
84  : getMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp")),
85  _dperm_or_cond_dvar(
86  _p_or_t == PorTchoice::pressure
87  ? getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")
88  : getMaterialProperty<std::vector<RealTensorValue>>(
89  "dPorousFlow_thermal_conductivity_qp_dvar"))
90 {
92  mooseError("PorousFlowPeacemanBorehole: You have specified function_of=porepressure, but you "
93  "do not have a quadpoint permeability material");
95  mooseError("PorousFlowPeacemanBorehole: You have specified function_of=temperature, but you do "
96  "not have a quadpoint thermal_conductivity material");
97 }
98 
99 void
101 {
103 
104  if (!_point_file.empty() && _zs[0] < _zs.back())
105  mooseError("PorousFlowPeacemanBorehole: The last entry in the point_file needs to be at the "
106  "bottom of the well_bore because this is the point where the function bottom_p_or_t "
107  "is evaluated. The depth of the first point is z=",
108  _zs[0],
109  " and the last point is z=",
110  _zs.back());
111 
112  // construct the rotation matrix needed to rotate the permeability
113  const unsigned int num_pts = _zs.size();
114  _rot_matrix.resize(std::max(num_pts - 1, (unsigned)1));
115  for (unsigned int i = 0; i + 1 < num_pts; ++i)
116  {
117  const RealVectorValue v2(_xs[i + 1] - _xs[i], _ys[i + 1] - _ys[i], _zs[i + 1] - _zs[i]);
119  }
120  if (num_pts == (unsigned)1)
122 }
123 
124 Real
126  const RealTensorValue & rot,
127  const Real & half_len,
128  const Elem * ele,
129  const Real & rad) const
130 // Peaceman's form for the borehole well constant
131 {
132  if (_well_constant > 0)
133  return _well_constant;
134 
135  // rot_perm has its "2" component lying along the half segment.
136  // We want to determine the eigenvectors of rot(0:1, 0:1), since, when
137  // rotated back to the original frame we will determine the element
138  // lengths along these directions
139  const RealTensorValue rot_perm = (rot * perm) * rot.transpose();
140  const Real trace2D = rot_perm(0, 0) + rot_perm(1, 1);
141  const Real det2D = rot_perm(0, 0) * rot_perm(1, 1) - rot_perm(0, 1) * rot_perm(1, 0);
142  const Real sq = std::sqrt(std::max(0.25 * trace2D * trace2D - det2D,
143  0.0)); // the std::max accounts for wierdo precision loss
144  const Real eig_val1 = 0.5 * trace2D + sq;
145  const Real eig_val2 = 0.5 * trace2D - sq;
146  RealVectorValue eig_vec1, eig_vec2;
147  if (sq > std::abs(trace2D) * 1E-7) // matrix is not a multiple of the identity (1E-7 accounts for
148  // precision in a crude way)
149  {
150  if (rot_perm(1, 0) != 0)
151  {
152  eig_vec1(0) = eig_val1 - rot_perm(1, 1);
153  eig_vec1(1) = rot_perm(1, 0);
154  eig_vec2(0) = eig_val2 - rot_perm(1, 1);
155  eig_vec2(1) = rot_perm(1, 0);
156  }
157  else if (rot_perm(0, 1) != 0)
158  {
159  eig_vec1(0) = rot_perm(0, 1);
160  eig_vec1(1) = eig_val1 - rot_perm(0, 0);
161  eig_vec2(0) = rot_perm(0, 1);
162  eig_vec2(1) = eig_val2 - rot_perm(0, 0);
163  }
164  else // off diagonal terms are both zero
165  {
166  eig_vec1(0) = 1.0;
167  eig_vec2(1) = 1.0;
168  }
169  }
170  else // matrix is basically a multiple of the identity
171  {
172  eig_vec1(0) = 1.0;
173  eig_vec2(1) = 1.0;
174  }
175 
176  // finally, rotate these to original frame and normalise
177  eig_vec1 = rot.transpose() * eig_vec1;
178  eig_vec1 /= std::sqrt(eig_vec1 * eig_vec1);
179  eig_vec2 = rot.transpose() * eig_vec2;
180  eig_vec2 /= std::sqrt(eig_vec2 * eig_vec2);
181 
182  // find the "length" of the element in these directions
183  // TODO - maybe better to use variance than max&min
184  Real max1 = eig_vec1 * ele->point(0);
185  Real max2 = eig_vec2 * ele->point(0);
186  Real min1 = max1;
187  Real min2 = max2;
188  Real proj;
189  for (unsigned int i = 1; i < ele->n_nodes(); i++)
190  {
191  proj = eig_vec1 * ele->point(i);
192  max1 = (max1 < proj) ? proj : max1;
193  min1 = (min1 < proj) ? min1 : proj;
194 
195  proj = eig_vec2 * ele->point(i);
196  max2 = (max2 < proj) ? proj : max2;
197  min2 = (min2 < proj) ? min2 : proj;
198  }
199  const Real ll1 = max1 - min1;
200  const Real ll2 = max2 - min2;
201 
202  Real r0;
203  if (eig_val1 <= 0.0)
204  r0 = _re_constant * ll1;
205  else if (eig_val2 <= 0.0)
206  r0 = _re_constant * ll2;
207  else
208  r0 = _re_constant *
209  std::sqrt(std::sqrt(eig_val1 / eig_val2) * std::pow(ll2, 2) +
210  std::sqrt(eig_val2 / eig_val1) * std::pow(ll1, 2)) /
211  (std::pow(eig_val1 / eig_val2, 0.25) + std::pow(eig_val2 / eig_val1, 0.25));
212 
213  const Real effective_perm = (det2D >= 0.0 ? std::sqrt(det2D) : 0.0);
214 
215  const Real halfPi = acos(0.0);
216 
217  if (r0 <= rad)
218  mooseError("The effective element size (about 0.2-times-true-ele-size) for an element "
219  "containing a Peaceman-type borehole must be (much) larger than the borehole radius "
220  "for the Peaceman formulation to be correct. Your element has effective size ",
221  r0,
222  " and the borehole radius is ",
223  rad,
224  "\n");
225 
226  return 4 * halfPi * effective_perm * half_len / std::log(r0 / rad);
227 }
228 
229 Real
230 PorousFlowPeacemanBorehole::computeQpBaseOutflow(unsigned current_dirac_ptid) const
231 {
232  const Real character = _character.value(_t, _q_point[_qp]);
233  if (character == 0.0)
234  return 0.0;
235 
236  const Real bh_pressure =
238  const Real pp = ptqp();
239 
240  Real outflow = 0.0; // this is the flow rate from porespace out of the system
241 
242  if (current_dirac_ptid > 0)
243  // contribution from half-segment "behind" this point (must have >1 point for
244  // current_dirac_ptid>0)
245  {
246  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
247  {
248  // injection, so outflow<0 || production, so outflow>0
249  const Real wc = wellConstant(_perm_or_cond[_qp],
250  _rot_matrix[current_dirac_ptid - 1],
251  _half_seg_len[current_dirac_ptid - 1],
253  _weight->at(current_dirac_ptid));
254  outflow += wc * (pp - bh_pressure);
255  }
256  }
257 
258  if (current_dirac_ptid + 1 < _zs.size() || _zs.size() == 1)
259  // contribution from half-segment "ahead of" this point, or we only have one point
260  {
261  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
262  {
263  // injection, so outflow<0 || // production, so outflow>0
264  const Real wc = wellConstant(_perm_or_cond[_qp],
265  _rot_matrix[current_dirac_ptid],
266  _half_seg_len[current_dirac_ptid],
268  _weight->at(current_dirac_ptid));
269  outflow += wc * (pp - bh_pressure);
270  }
271  }
272 
273  return outflow * _test[_i][_qp] * std::abs(character);
274 }
275 
276 void
278  unsigned current_dirac_ptid,
279  Real & outflow,
280  Real & outflowp) const
281 {
282  outflow = 0.0;
283  outflowp = 0.0;
284 
285  const Real character = _character.value(_t, _q_point[_qp]);
286  if (character == 0.0)
287  return;
288 
290  return;
291  const unsigned pvar = _dictator.porousFlowVariableNum(jvar);
292 
293  const Real bh_pressure =
295  const Real pp = ptqp();
296  const Real pp_prime = dptqp(pvar) * _phi[_j][_qp];
297 
298  if (current_dirac_ptid > 0)
299  // contribution from half-segment "behind" this point
300  {
301  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
302  {
303  // injection, so outflow<0 || // production, so outflow>0
304  const Real wc = wellConstant(_perm_or_cond[_qp],
305  _rot_matrix[current_dirac_ptid - 1],
306  _half_seg_len[current_dirac_ptid - 1],
308  _weight->at(current_dirac_ptid));
309  outflowp += wc * pp_prime;
310  outflow += wc * (pp - bh_pressure);
311  }
312  }
313 
314  if (current_dirac_ptid < _zs.size() - 1 || _zs.size() == 1)
315  // contribution from half-segment "ahead of" this point
316  {
317  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
318  {
319  // injection, so outflow<0 || // production, so outflow>0
320  const Real wc = wellConstant(_perm_or_cond[_qp],
321  _rot_matrix[current_dirac_ptid],
322  _half_seg_len[current_dirac_ptid],
324  _weight->at(current_dirac_ptid));
325  outflowp += wc * pp_prime;
326  outflow += wc * (pp - bh_pressure);
327  }
328  }
329 
330  outflowp *= _test[_i][_qp] * std::abs(character);
331  outflow *= _test[_i][_qp] * std::abs(character);
332 }
PorousFlowPeacemanBorehole(const InputParameters &parameters)
const std::vector< Real > *const _weight
void computeQpBaseOutflowJacobian(unsigned jvar, unsigned current_dirac_ptid, Real &outflow, Real &outflowp) const override
Calculates the BaseOutflow as well as its derivative wrt jvar. Derived classes should override this...
bool notPorousFlowVariable(unsigned int moose_var_num) const
Returns true if moose_var_num is not a porous flow variabe.
const MooseArray< Point > & _q_point
const PorousFlowDictator & _dictator
PorousFlowDictator UserObject.
static InputParameters validParams()
Real computeQpBaseOutflow(unsigned current_dirac_ptid) const override
Returns the flux from the line sink (before modification by mobility, etc). Derived classes should ov...
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
void mooseError(Args &&... args)
const Function & _character
If positive then the borehole acts as a sink (producion well) for porepressure > borehole pressure...
const bool _has_permeability
Whether there is a quadpoint permeability material (for error checking)
registerMooseObject("PorousFlowApp", PorousFlowPeacemanBorehole)
const OutputTools< T >::VariableTestValue & _test
const OutputTools< T >::VariablePhiValue & _phi
std::vector< Real > _ys
y points of the borehole
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
const Function & _p_bot
Bottomhole pressure of borehole.
void addRequiredParam(const std::string &name, const std::string &doc_string)
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
TensorValue< Real > RealTensorValue
std::vector< Real > _half_seg_len
0.5*(length of polyline segments between points)
const RealVectorValue _line_direction
Line direction. This is only used if there is only one borehole point.
Real wellConstant(const RealTensorValue &perm, const RealTensorValue &rot, const Real &half_len, const Elem *ele, const Real &rad) const
Calculates Peaceman&#39;s form of the borehole well constant Z Chen, Y Zhang, Well flow models for variou...
virtual void initialSetup() override
Real ptqp() const
If _p_or_t==0, then returns the quadpoint porepressure, else returns the quadpoint temperature...
std::vector< Real > _zs
z points of borehole
Approximates a borehole by a sequence of Dirac Points.
std::vector< RealTensorValue > _rot_matrix
Rotation matrix used in well_constant calculation.
static InputParameters validParams()
Creates a new PorousFlowPeacemanBorehole This reads the file containing the lines of the form radius ...
const Elem *const & _current_elem
PorTchoice
whether the flux is a function of pressure or temperature
GenericRealTensorValue< is_ad > rotVecToZ(GenericRealVectorValue< is_ad > vec)
const std::string _point_file
File defining the geometry of the borehole.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int _j
unsigned int _qp
static const std::string pressure
Definition: NS.h:56
void mooseError(Args &&... args) const
virtual void initialSetup() override
const MaterialProperty< RealTensorValue > & _perm_or_cond
Permeability or conductivity of porous material.
void addClassDescription(const std::string &doc_string)
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
The PorousFlow variable number.
Approximates a line sink a sequence of Dirac Points.
std::vector< Real > _xs
x points of the borehole
virtual Real value(Real t, const Point &p) const
Point _bottom_point
The bottom point of the borehole (where bottom_pressure is defined)
unsigned int _i
const bool _has_thermal_conductivity
Whether there is a quadpoint thermal conductivity material (for error checking)
enum PorousFlowLineSink::PorTchoice _p_or_t
MooseUnits pow(const MooseUnits &, int)
Real dptqp(unsigned pvar) const
If _p_or_t==0, then returns d(quadpoint porepressure)/d(PorousFlow variable), else returns d(quadpoin...
const RealVectorValue _unit_weight
Unit weight of fluid in borehole (for calculating bottomhole pressure at each Dirac Point) ...