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 " 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 " 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");
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.)");
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");
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 " 70 _character(getFunction(
"character")),
71 _p_bot(getFunction(
"bottom_p_or_t")),
73 _re_constant(getParam<
Real>(
"re_constant")),
74 _well_constant(getParam<
Real>(
"well_constant")),
77 hasMaterialProperty<
std::vector<
RealTensorValue>>(
"dPorousFlow_permeability_qp_dvar")),
78 _has_thermal_conductivity(
79 hasMaterialProperty<
RealTensorValue>(
"PorousFlow_thermal_conductivity_qp") &&
81 "dPorousFlow_thermal_conductivity_qp_dvar")),
84 : getMaterialProperty<
RealTensorValue>(
"PorousFlow_thermal_conductivity_qp")),
87 ? getMaterialProperty<
std::vector<
RealTensorValue>>(
"dPorousFlow_permeability_qp_dvar")
89 "dPorousFlow_thermal_conductivity_qp_dvar"))
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");
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=",
109 " and the last point is z=",
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)
120 if (num_pts == (
unsigned)1)
127 const Real & half_len,
129 const Real & rad)
const 132 if (_well_constant > 0)
133 return _well_constant;
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,
144 const Real eig_val1 = 0.5 * trace2D + sq;
145 const Real eig_val2 = 0.5 * trace2D - sq;
150 if (rot_perm(1, 0) != 0)
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);
157 else if (rot_perm(0, 1) != 0)
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);
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);
184 Real max1 = eig_vec1 * ele->point(0);
185 Real max2 = eig_vec2 * ele->point(0);
189 for (
unsigned int i = 1; i < ele->n_nodes(); i++)
191 proj = eig_vec1 * ele->point(i);
192 max1 = (max1 < proj) ? proj : max1;
193 min1 = (min1 < proj) ? min1 : proj;
195 proj = eig_vec2 * ele->point(i);
196 max2 = (max2 < proj) ? proj : max2;
197 min2 = (min2 < proj) ? min2 : proj;
199 const Real ll1 = max1 - min1;
200 const Real ll2 = max2 - min2;
204 r0 = _re_constant * ll1;
205 else if (eig_val2 <= 0.0)
206 r0 = _re_constant * ll2;
211 (
std::pow(eig_val1 / eig_val2, 0.25) +
std::pow(eig_val2 / eig_val1, 0.25));
213 const Real effective_perm = (det2D >= 0.0 ?
std::sqrt(det2D) : 0.0);
215 const Real halfPi = acos(0.0);
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 ",
222 " and the borehole radius is ",
226 return 4 * halfPi * effective_perm * half_len / std::log(r0 / rad);
233 if (character == 0.0)
236 const Real bh_pressure =
242 if (current_dirac_ptid > 0)
246 if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
253 _weight->at(current_dirac_ptid));
254 outflow += wc * (pp - bh_pressure);
258 if (current_dirac_ptid + 1 <
_zs.size() ||
_zs.size() == 1)
261 if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
268 _weight->at(current_dirac_ptid));
269 outflow += wc * (pp - bh_pressure);
278 unsigned current_dirac_ptid,
280 Real & outflowp)
const 286 if (character == 0.0)
293 const Real bh_pressure =
298 if (current_dirac_ptid > 0)
301 if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
308 _weight->at(current_dirac_ptid));
309 outflowp += wc * pp_prime;
310 outflow += wc * (pp - bh_pressure);
314 if (current_dirac_ptid <
_zs.size() - 1 ||
_zs.size() == 1)
317 if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
324 _weight->at(current_dirac_ptid));
325 outflowp += wc * pp_prime;
326 outflow += wc * (pp - bh_pressure);
PorousFlowPeacemanBorehole(const InputParameters ¶meters)
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 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.
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'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
static const std::string pressure
void mooseError(Args &&... args) const
virtual void initialSetup() override
const MaterialProperty< RealTensorValue > & _perm_or_cond
Permeability or conductivity of porous material.
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)
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) ...