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 "PorousFlowPeacemanBorehole.h"
11 : #include "RotationMatrix.h"
12 : #include "Function.h"
13 :
14 : registerMooseObject("PorousFlowApp", PorousFlowPeacemanBorehole);
15 :
16 : InputParameters
17 1076 : PorousFlowPeacemanBorehole::validParams()
18 : {
19 1076 : InputParameters params = PorousFlowLineSink::validParams();
20 2152 : 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 2152 : 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.");
34 2152 : params.addRequiredParam<RealVectorValue>(
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 2152 : params.addParam<Real>("re_constant",
44 2152 : 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 2152 : params.addParam<Real>("well_constant",
52 2152 : -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 1076 : 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 1076 : return params;
66 0 : }
67 :
68 610 : PorousFlowPeacemanBorehole::PorousFlowPeacemanBorehole(const InputParameters & parameters)
69 : : PorousFlowLineSink(parameters),
70 588 : _character(getFunction("character")),
71 588 : _p_bot(getFunction("bottom_p_or_t")),
72 1176 : _unit_weight(getParam<RealVectorValue>("unit_weight")),
73 1176 : _re_constant(getParam<Real>("re_constant")),
74 1176 : _well_constant(getParam<Real>("well_constant")),
75 588 : _has_permeability(
76 588 : hasMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp") &&
77 1174 : hasMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
78 588 : _has_thermal_conductivity(
79 588 : hasMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp") &&
80 984 : hasMaterialProperty<std::vector<RealTensorValue>>(
81 : "dPorousFlow_thermal_conductivity_qp_dvar")),
82 1176 : _perm_or_cond(_p_or_t == PorTchoice::pressure
83 588 : ? getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")
84 724 : : getMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp")),
85 588 : _dperm_or_cond_dvar(
86 588 : _p_or_t == PorTchoice::pressure
87 1176 : ? getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")
88 724 : : getMaterialProperty<std::vector<RealTensorValue>>(
89 610 : "dPorousFlow_thermal_conductivity_qp_dvar"))
90 : {
91 588 : if (_p_or_t == PorTchoice::pressure && !_has_permeability)
92 2 : mooseError("PorousFlowPeacemanBorehole: You have specified function_of=porepressure, but you "
93 : "do not have a quadpoint permeability material");
94 586 : if (_p_or_t == PorTchoice::temperature && !_has_thermal_conductivity)
95 2 : mooseError("PorousFlowPeacemanBorehole: You have specified function_of=temperature, but you do "
96 : "not have a quadpoint thermal_conductivity material");
97 584 : }
98 :
99 : void
100 581 : PorousFlowPeacemanBorehole::initialSetup()
101 : {
102 581 : PorousFlowLineGeometry::initialSetup();
103 :
104 577 : if (!_point_file.empty() && _zs[0] < _zs.back())
105 2 : 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 575 : const unsigned int num_pts = _zs.size();
114 839 : _rot_matrix.resize(std::max(num_pts - 1, (unsigned)1));
115 1062 : for (unsigned int i = 0; i + 1 < num_pts; ++i)
116 : {
117 487 : const RealVectorValue v2(_xs[i + 1] - _xs[i], _ys[i + 1] - _ys[i], _zs[i + 1] - _zs[i]);
118 487 : _rot_matrix[i] = RotationMatrix::rotVecToZ(v2);
119 : }
120 575 : if (num_pts == (unsigned)1)
121 264 : _rot_matrix[0] = RotationMatrix::rotVecToZ(_line_direction);
122 575 : }
123 :
124 : Real
125 1012258 : PorousFlowPeacemanBorehole::wellConstant(const RealTensorValue & perm,
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 1012258 : 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 1012258 : const RealTensorValue rot_perm = (rot * perm) * rot.transpose();
140 1012258 : const Real trace2D = rot_perm(0, 0) + rot_perm(1, 1);
141 1012258 : const Real det2D = rot_perm(0, 0) * rot_perm(1, 1) - rot_perm(0, 1) * rot_perm(1, 0);
142 2024516 : const Real sq = std::sqrt(std::max(0.25 * trace2D * trace2D - det2D,
143 1012258 : 0.0)); // the std::max accounts for wierdo precision loss
144 1012258 : const Real eig_val1 = 0.5 * trace2D + sq;
145 1012258 : const Real eig_val2 = 0.5 * trace2D - sq;
146 : RealVectorValue eig_vec1, eig_vec2;
147 1012258 : 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 148360 : if (rot_perm(1, 0) != 0)
151 : {
152 120640 : eig_vec1(0) = eig_val1 - rot_perm(1, 1);
153 120640 : eig_vec1(1) = rot_perm(1, 0);
154 120640 : eig_vec2(0) = eig_val2 - rot_perm(1, 1);
155 120640 : eig_vec2(1) = rot_perm(1, 0);
156 : }
157 27720 : else if (rot_perm(0, 1) != 0)
158 : {
159 0 : eig_vec1(0) = rot_perm(0, 1);
160 0 : eig_vec1(1) = eig_val1 - rot_perm(0, 0);
161 0 : eig_vec2(0) = rot_perm(0, 1);
162 0 : eig_vec2(1) = eig_val2 - rot_perm(0, 0);
163 : }
164 : else // off diagonal terms are both zero
165 : {
166 27720 : eig_vec1(0) = 1.0;
167 27720 : eig_vec2(1) = 1.0;
168 : }
169 : }
170 : else // matrix is basically a multiple of the identity
171 : {
172 863898 : eig_vec1(0) = 1.0;
173 863898 : eig_vec2(1) = 1.0;
174 : }
175 :
176 : // finally, rotate these to original frame and normalise
177 1012258 : eig_vec1 = rot.transpose() * eig_vec1;
178 1012258 : eig_vec1 /= std::sqrt(eig_vec1 * eig_vec1);
179 1012258 : eig_vec2 = rot.transpose() * eig_vec2;
180 1012258 : 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 8098064 : for (unsigned int i = 1; i < ele->n_nodes(); i++)
190 : {
191 : proj = eig_vec1 * ele->point(i);
192 7085806 : max1 = (max1 < proj) ? proj : max1;
193 7085806 : min1 = (min1 < proj) ? min1 : proj;
194 :
195 : proj = eig_vec2 * ele->point(i);
196 7085806 : max2 = (max2 < proj) ? proj : max2;
197 7085806 : min2 = (min2 < proj) ? min2 : proj;
198 : }
199 1012258 : const Real ll1 = max1 - min1;
200 1012258 : const Real ll2 = max2 - min2;
201 :
202 : Real r0;
203 1012258 : if (eig_val1 <= 0.0)
204 0 : r0 = _re_constant * ll1;
205 1012258 : else if (eig_val2 <= 0.0)
206 18560 : r0 = _re_constant * ll2;
207 : else
208 993698 : r0 = _re_constant *
209 993698 : std::sqrt(std::sqrt(eig_val1 / eig_val2) * std::pow(ll2, 2) +
210 993698 : std::sqrt(eig_val2 / eig_val1) * std::pow(ll1, 2)) /
211 993698 : (std::pow(eig_val1 / eig_val2, 0.25) + std::pow(eig_val2 / eig_val1, 0.25));
212 :
213 1012258 : const Real effective_perm = (det2D >= 0.0 ? std::sqrt(det2D) : 0.0);
214 :
215 : const Real halfPi = acos(0.0);
216 :
217 1012258 : if (r0 <= rad)
218 2 : 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 1012256 : return 4 * halfPi * effective_perm * half_len / std::log(r0 / rad);
227 : }
228 :
229 : Real
230 222018 : PorousFlowPeacemanBorehole::computeQpBaseOutflow(unsigned current_dirac_ptid) const
231 : {
232 222018 : const Real character = _character.value(_t, _q_point[_qp]);
233 222018 : if (character == 0.0)
234 : return 0.0;
235 :
236 : const Real bh_pressure =
237 186018 : _p_bot.value(_t, _bottom_point) + _unit_weight * (_q_point[_qp] - _bottom_point);
238 186018 : const Real pp = ptqp();
239 :
240 : Real outflow = 0.0; // this is the flow rate from porespace out of the system
241 :
242 186018 : 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 77288 : if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
247 : {
248 : // injection, so outflow<0 || production, so outflow>0
249 77144 : const Real wc = wellConstant(_perm_or_cond[_qp],
250 : _rot_matrix[current_dirac_ptid - 1],
251 77144 : _half_seg_len[current_dirac_ptid - 1],
252 77144 : _current_elem,
253 77144 : _weight->at(current_dirac_ptid));
254 77144 : outflow += wc * (pp - bh_pressure);
255 : }
256 : }
257 :
258 186018 : 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 123130 : if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
262 : {
263 : // injection, so outflow<0 || // production, so outflow>0
264 122506 : const Real wc = wellConstant(_perm_or_cond[_qp],
265 : _rot_matrix[current_dirac_ptid],
266 : _half_seg_len[current_dirac_ptid],
267 122506 : _current_elem,
268 122506 : _weight->at(current_dirac_ptid));
269 122504 : outflow += wc * (pp - bh_pressure);
270 : }
271 : }
272 :
273 186016 : return outflow * _test[_i][_qp] * std::abs(character);
274 : }
275 :
276 : void
277 838456 : PorousFlowPeacemanBorehole::computeQpBaseOutflowJacobian(unsigned jvar,
278 : unsigned current_dirac_ptid,
279 : Real & outflow,
280 : Real & outflowp) const
281 : {
282 838456 : outflow = 0.0;
283 838456 : outflowp = 0.0;
284 :
285 838456 : const Real character = _character.value(_t, _q_point[_qp]);
286 838456 : if (character == 0.0)
287 : return;
288 :
289 804856 : if (_dictator.notPorousFlowVariable(jvar))
290 : return;
291 804856 : const unsigned pvar = _dictator.porousFlowVariableNum(jvar);
292 :
293 : const Real bh_pressure =
294 804856 : _p_bot.value(_t, _bottom_point) + _unit_weight * (_q_point[_qp] - _bottom_point);
295 804856 : const Real pp = ptqp();
296 804856 : const Real pp_prime = dptqp(pvar) * _phi[_j][_qp];
297 :
298 804856 : if (current_dirac_ptid > 0)
299 : // contribution from half-segment "behind" this point
300 : {
301 387328 : if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
302 : {
303 : // injection, so outflow<0 || // production, so outflow>0
304 385024 : const Real wc = wellConstant(_perm_or_cond[_qp],
305 : _rot_matrix[current_dirac_ptid - 1],
306 385024 : _half_seg_len[current_dirac_ptid - 1],
307 385024 : _current_elem,
308 385024 : _weight->at(current_dirac_ptid));
309 385024 : outflowp += wc * pp_prime;
310 385024 : outflow += wc * (pp - bh_pressure);
311 : }
312 : }
313 :
314 804856 : if (current_dirac_ptid < _zs.size() - 1 || _zs.size() == 1)
315 : // contribution from half-segment "ahead of" this point
316 : {
317 430968 : if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
318 : {
319 : // injection, so outflow<0 || // production, so outflow>0
320 427584 : const Real wc = wellConstant(_perm_or_cond[_qp],
321 : _rot_matrix[current_dirac_ptid],
322 : _half_seg_len[current_dirac_ptid],
323 427584 : _current_elem,
324 427584 : _weight->at(current_dirac_ptid));
325 427584 : outflowp += wc * pp_prime;
326 427584 : outflow += wc * (pp - bh_pressure);
327 : }
328 : }
329 :
330 804856 : outflowp *= _test[_i][_qp] * std::abs(character);
331 804856 : outflow *= _test[_i][_qp] * std::abs(character);
332 : }
|
