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