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