www.mooseframework.org
PeacemanBorehole.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 "PeacemanBorehole.h"
9 #include "RotationMatrix.h"
10 
11 #include <fstream>
12 
13 template <>
14 InputParameters
16 {
17  InputParameters params = validParams<DiracKernel>();
18  params.addRequiredParam<FunctionName>(
19  "character",
20  "If zero then borehole does nothing. If positive the borehole acts as a sink "
21  "(production well) for porepressure > borehole pressure, and does nothing "
22  "otherwise. If negative the borehole acts as a source (injection well) for "
23  "porepressure < borehole pressure, and does nothing otherwise. The flow rate "
24  "to/from the borehole is multiplied by |character|, so usually character = +/- "
25  "1, but you can specify other quantities to provide an overall scaling to the "
26  "flow if you like.");
27  params.addRequiredParam<Real>("bottom_pressure", "Pressure at the bottom of the borehole");
28  params.addRequiredParam<RealVectorValue>(
29  "unit_weight",
30  "(fluid_density*gravitational_acceleration) as a vector pointing downwards. "
31  "Note that the borehole pressure at a given z position is bottom_pressure + "
32  "unit_weight*(p - p_bottom), where p=(x,y,z) and p_bottom=(x,y,z) of the "
33  "bottom point of the borehole. If you don't want bottomhole pressure to vary "
34  "in the borehole just set unit_weight=0. Typical value is = (0,0,-1E4)");
35  params.addRequiredParam<std::string>(
36  "point_file",
37  "The file containing the borehole radii and coordinates of the point sinks "
38  "that approximate the borehole. Each line in the file must contain a "
39  "space-separated radius and coordinate. Ie r x y z. The last point in the "
40  "file is defined as the borehole bottom, where the borehole pressure is "
41  "bottom_pressure. If your file contains just one point, you must also specify "
42  "the borehole_length and borehole_direction. Note that you will get "
43  "segementation faults if your points do not lie within your mesh!");
44  params.addRequiredParam<UserObjectName>(
45  "SumQuantityUO",
46  "User Object of type=RichardsSumQuantity in which to place the total "
47  "outflow from the borehole for each time step.");
48  params.addParam<Real>("re_constant",
49  0.28,
50  "The dimensionless constant used in evaluating the borehole effective "
51  "radius. This depends on the meshing scheme. Peacemann "
52  "finite-difference calculations give 0.28, while for rectangular finite "
53  "elements the result is closer to 0.1594. (See Eqn(4.13) of Z Chen, Y "
54  "Zhang, Well flow models for various numerical methods, Int J Num "
55  "Analysis and Modeling, 3 (2008) 375-388.)");
56  params.addParam<Real>("well_constant",
57  -1.0,
58  "Usually this is calculated internally from the element geometry, the "
59  "local borehole direction and segment length, and the permeability. "
60  "However, if this parameter is given as a positive number then this "
61  "number is used instead of the internal calculation. This speeds up "
62  "computation marginally. re_constant becomes irrelevant");
63  params.addRangeCheckedParam<Real>(
64  "borehole_length",
65  0.0,
66  "borehole_length>=0",
67  "Borehole length. Note this is only used if there is only one point in the point_file.");
68  params.addParam<RealVectorValue>(
69  "borehole_direction",
70  RealVectorValue(0, 0, 1),
71  "Borehole direction. Note this is only used if there is only one point in the point_file.");
72  params.addClassDescription("Approximates a borehole in the mesh using the Peaceman approach, ie "
73  "using a number of point sinks with given radii whose positions are "
74  "read from a file");
75  return params;
76 }
77 
78 PeacemanBorehole::PeacemanBorehole(const InputParameters & parameters)
79  : DiracKernel(parameters),
80  _re_constant(getParam<Real>("re_constant")),
81  _well_constant(getParam<Real>("well_constant")),
82  _borehole_length(getParam<Real>("borehole_length")),
83  _borehole_direction(getParam<RealVectorValue>("borehole_direction")),
84  _point_file(getParam<std::string>("point_file")),
85  _character(getFunction("character")),
86  _p_bot(getParam<Real>("bottom_pressure")),
87  _unit_weight(getParam<RealVectorValue>("unit_weight")),
88  _total_outflow_mass(
89  const_cast<RichardsSumQuantity &>(getUserObject<RichardsSumQuantity>("SumQuantityUO")))
90 {
91  // zero the outflow mass
93 
94  // open file
95  std::ifstream file(_point_file.c_str());
96  if (!file.good())
97  mooseError("Error opening file '" + _point_file + "' from a Peaceman-type Borehole.");
98 
99  // construct the arrays of radius, x, y and z
100  std::vector<Real> scratch;
101  while (parseNextLineReals(file, scratch))
102  {
103  if (scratch.size() >= 2)
104  {
105  _rs.push_back(scratch[0]);
106  _xs.push_back(scratch[1]);
107  if (scratch.size() >= 3)
108  _ys.push_back(scratch[2]);
109  else
110  _ys.push_back(0.0);
111  if (scratch.size() >= 4)
112  _zs.push_back(scratch[3]);
113  else
114  _zs.push_back(0.0);
115  }
116  }
117 
118  file.close();
119 
120  const int num_pts = _zs.size();
121  _bottom_point(0) = _xs[num_pts - 1];
122  _bottom_point(1) = _ys[num_pts - 1];
123  _bottom_point(2) = _zs[num_pts - 1];
124 
125  // construct the line-segment lengths between each point
126  _half_seg_len.resize(std::max(num_pts - 1, 1));
127  for (unsigned int i = 0; i + 1 < _xs.size(); ++i)
128  {
129  _half_seg_len[i] =
130  0.5 * std::sqrt(std::pow(_xs[i + 1] - _xs[i], 2) + std::pow(_ys[i + 1] - _ys[i], 2) +
131  std::pow(_zs[i + 1] - _zs[i], 2));
132  if (_half_seg_len[i] == 0)
133  mooseError("Peaceman-type borehole has a zero-segment length at (x,y,z) = ",
134  _xs[i],
135  " ",
136  _ys[i],
137  " ",
138  _zs[i],
139  "\n");
140  }
141  if (num_pts == 1)
143 
144  // construct the rotation matrix needed to rotate the permeability
145  _rot_matrix.resize(std::max(num_pts - 1, 1));
146  for (unsigned int i = 0; i + 1 < _xs.size(); ++i)
147  {
148  const RealVectorValue v2(_xs[i + 1] - _xs[i], _ys[i + 1] - _ys[i], _zs[i + 1] - _zs[i]);
149  _rot_matrix[i] = RotationMatrix::rotVecToZ(v2);
150  }
151  if (num_pts == 1)
152  _rot_matrix[0] = RotationMatrix::rotVecToZ(_borehole_direction);
153 }
154 
155 bool
156 PeacemanBorehole::parseNextLineReals(std::ifstream & ifs, std::vector<Real> & myvec)
157 // reads a space-separated line of floats from ifs and puts in myvec
158 {
159  std::string line;
160  myvec.clear();
161  bool gotline(false);
162  if (getline(ifs, line))
163  {
164  gotline = true;
165 
166  // Harvest floats separated by whitespace
167  std::istringstream iss(line);
168  Real f;
169  while (iss >> f)
170  {
171  myvec.push_back(f);
172  }
173  }
174  return gotline;
175 }
176 
177 void
179 {
180  // This function gets called just before the DiracKernel is evaluated
181  // so this is a handy place to zero this out.
183 
184  // Add point using the unique ID "i", let the DiracKernel take
185  // care of the caching. This should be fast after the first call,
186  // as long as the points don't move around.
187  for (unsigned int i = 0; i < _zs.size(); i++)
188  addPoint(Point(_xs[i], _ys[i], _zs[i]), i);
189 }
190 
191 Real
192 PeacemanBorehole::wellConstant(const RealTensorValue & perm,
193  const RealTensorValue & rot,
194  const Real & half_len,
195  const Elem * ele,
196  const Real & rad)
197 // Peaceman's form for the borehole well constant
198 {
199  if (_well_constant > 0)
200  return _well_constant;
201 
202  // rot_perm has its "2" component lying along the half segment
203  // we want to determine the eigenvectors of rot(0:1, 0:1), since, when
204  // rotated back to the original frame we will determine the element
205  // lengths along these directions
206  const RealTensorValue rot_perm = (rot * perm) * rot.transpose();
207  const Real trace2D = rot_perm(0, 0) + rot_perm(1, 1);
208  const Real det2D = rot_perm(0, 0) * rot_perm(1, 1) - rot_perm(0, 1) * rot_perm(1, 0);
209  const Real sq = std::sqrt(std::max(0.25 * trace2D * trace2D - det2D,
210  0.0)); // the std::max accounts for wierdo precision loss
211  const Real eig_val1 = 0.5 * trace2D + sq;
212  const Real eig_val2 = 0.5 * trace2D - sq;
213  RealVectorValue eig_vec1, eig_vec2;
214  if (sq > std::abs(trace2D) * 1E-7) // matrix is not a multiple of the identity (1E-7 accounts for
215  // precision in a crude way)
216  {
217  if (rot_perm(1, 0) != 0)
218  {
219  eig_vec1(0) = eig_val1 - rot_perm(1, 1);
220  eig_vec1(1) = rot_perm(1, 0);
221  eig_vec2(0) = eig_val2 - rot_perm(1, 1);
222  eig_vec2(1) = rot_perm(1, 0);
223  }
224  else if (rot_perm(0, 1) != 0)
225  {
226  eig_vec1(0) = rot_perm(0, 1);
227  eig_vec1(1) = eig_val1 - rot_perm(0, 0);
228  eig_vec2(0) = rot_perm(0, 1);
229  eig_vec2(1) = eig_val2 - rot_perm(0, 0);
230  }
231  else // off diagonal terms are both zero
232  {
233  eig_vec1(0) = 1;
234  eig_vec2(1) = 1;
235  }
236  }
237  else // matrix is basically a multiple of the identity
238  {
239  eig_vec1(0) = 1;
240  eig_vec2(1) = 1;
241  }
242 
243  // finally, rotate these to original frame and normalise
244  eig_vec1 = rot.transpose() * eig_vec1;
245  eig_vec1 /= std::sqrt(eig_vec1 * eig_vec1);
246  eig_vec2 = rot.transpose() * eig_vec2;
247  eig_vec2 /= std::sqrt(eig_vec2 * eig_vec2);
248 
249  // find the "length" of the element in these directions
250  // TODO - probably better to use variance than max&min
251  Real max1 = eig_vec1 * ele->point(0);
252  Real max2 = eig_vec2 * ele->point(0);
253  Real min1 = max1;
254  Real min2 = max2;
255  Real proj;
256  for (unsigned int i = 1; i < ele->n_nodes(); i++)
257  {
258  proj = eig_vec1 * ele->point(i);
259  max1 = (max1 < proj) ? proj : max1;
260  min1 = (min1 < proj) ? min1 : proj;
261 
262  proj = eig_vec2 * ele->point(i);
263  max2 = (max2 < proj) ? proj : max2;
264  min2 = (min2 < proj) ? min2 : proj;
265  }
266  const Real ll1 = max1 - min1;
267  const Real ll2 = max2 - min2;
268 
269  const Real r0 = _re_constant * std::sqrt(std::sqrt(eig_val1 / eig_val2) * std::pow(ll2, 2) +
270  std::sqrt(eig_val2 / eig_val1) * std::pow(ll1, 2)) /
271  (std::pow(eig_val1 / eig_val2, 0.25) + std::pow(eig_val2 / eig_val1, 0.25));
272 
273  const Real effective_perm = std::sqrt(det2D);
274 
275  const Real halfPi = acos(0.0);
276 
277  if (r0 <= rad)
278  mooseError("The effective element size (about 0.2-times-true-ele-size) for an element "
279  "containing a Peaceman-type borehole must be (much) larger than the borehole radius "
280  "for the Peaceman formulation to be correct. Your element has effective size ",
281  r0,
282  " and the borehole radius is ",
283  rad,
284  "\n");
285 
286  return 4 * halfPi * effective_perm * half_len / std::log(r0 / rad);
287 }
const Real _well_constant
well constant
const RealVectorValue _borehole_direction
borehole direction. Note this is only used if there is only one borehole point
Sums into _total This is used, for instance, to record the total mass flowing into a borehole...
const Real _re_constant
borehole constant
PeacemanBorehole(const InputParameters &parameters)
Creates a new PeacemanBorehole This reads the file containing the lines of the form radius x y z that...
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...
InputParameters validParams< PeacemanBorehole >()
std::vector< Real > _rs
radii of the borehole
bool parseNextLineReals(std::ifstream &ifs, std::vector< Real > &myvec)
reads a space-separated line of floats from ifs and puts in myvec
std::vector< Real > _xs
x points of the borehole
std::vector< Real > _ys
y points of the borehole
void zero()
sets _total = 0
RichardsSumQuantity & _total_outflow_mass
This is used to hold the total fluid flowing into the borehole Hence, it is positive for production w...
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
const Real _borehole_length
borehole length. Note this is only used if there is only one borehole point
virtual void addPoints()
Add Dirac Points to the borehole.
const std::string _point_file
File defining the geometry of the borehole.
Point _bottom_point
the bottom point of the borehole (where bottom_pressure is defined)
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