Approximates a borehole by a sequence of Dirac Points. More...

#include <PeacemanBorehole.h>

Inheritance diagram for PeacemanBorehole:

Public Member Functions
	PeacemanBorehole (const InputParameters &parameters)
	Creates a new PeacemanBorehole This reads the file containing the lines of the form radius x y z that defines the borehole geometry. More...

Protected Member Functions
virtual void	addPoints ()
	Add Dirac Points to the borehole. More...

bool	parseNextLineReals (std::ifstream &ifs, std::vector< Real > &myvec)
	reads a space-separated line of floats from ifs and puts in myvec More...

Real	wellConstant (const RealTensorValue &perm, const RealTensorValue &rot, const Real &half_len, const Elem *ele, const Real &rad)
	Calculates Peaceman's form of the borehole well constant Z Chen, Y Zhang, Well flow models for various numerical methods, Int J Num Analysis and Modeling, 3 (2008) 375-388. More...

Protected Attributes
const Function &	_character
	If positive then the borehole acts as a sink (producion well) for porepressure > borehole pressure, and does nothing otherwise If negative then the borehole acts as a source (injection well) for porepressure < borehole pressure, and does nothing otherwise The flow rate to/from the borehole is multiplied by \|character\|, so usually character = +/- 1. More...

const Real	_p_bot
	bottomhole pressure of borehole More...

const RealVectorValue	_unit_weight
	unit weight of fluid in borehole (for calculating bottomhole pressure at each Dirac Point) More...

RichardsSumQuantity &	_total_outflow_mass
	This is used to hold the total fluid flowing into the borehole Hence, it is positive for production wells where fluid is flowing from porespace into the borehole and removed from the model. More...

std::vector< Real >	_rs
	radii of the borehole More...

std::vector< Real >	_xs
	x points of the borehole More...

std::vector< Real >	_ys
	y points of the borehole More...

std::vector< Real >	_zs
	z points of borehole More...

Point	_bottom_point
	the bottom point of the borehole (where bottom_pressure is defined) More...

std::vector< Real >	_half_seg_len
	0.5*(length of polyline segments between points) More...

std::vector< RealTensorValue >	_rot_matrix
	rotation matrix used in well_constant calculation More...

Private Attributes
const Real	_re_constant
	borehole constant More...

const Real	_well_constant
	well constant More...

const Real	_borehole_length
	borehole length. Note this is only used if there is only one borehole point More...

const RealVectorValue	_borehole_direction
	borehole direction. Note this is only used if there is only one borehole point More...

const std::string	_point_file
	File defining the geometry of the borehole. More...

Detailed Description

Approximates a borehole by a sequence of Dirac Points.

Definition at line 25 of file PeacemanBorehole.h.

Constructor & Destructor Documentation

◆ PeacemanBorehole()

PeacemanBorehole::PeacemanBorehole ( const InputParameters & parameters )

Creates a new PeacemanBorehole This reads the file containing the lines of the form radius x y z that defines the borehole geometry.

It also calculates segment-lengths and rotation matrices needed for computing the borehole well constant

Definition at line 80 of file PeacemanBorehole.C.

   : DiracKernel(parameters),
     _re_constant(getParam<Real>("re_constant")),
     _well_constant(getParam<Real>("well_constant")),
     _borehole_length(getParam<Real>("borehole_length")),
     _borehole_direction(getParam<RealVectorValue>("borehole_direction")),
     _point_file(getParam<std::string>("point_file")),
     _character(getFunction("character")),
     _p_bot(getParam<Real>("bottom_pressure")),
     _unit_weight(getParam<RealVectorValue>("unit_weight")),
     _total_outflow_mass(
         const_cast<RichardsSumQuantity &>(getUserObject<RichardsSumQuantity>("SumQuantityUO")))
 {
   // zero the outflow mass
   _total_outflow_mass.zero();
  
   // open file
   std::ifstream file(_point_file.c_str());
   if (!file.good())
     mooseError("Error opening file '" + _point_file + "' from a Peaceman-type Borehole.");
  
   // construct the arrays of radius, x, y and z
   std::vector<Real> scratch;
   while (parseNextLineReals(file, scratch))
   {
     if (scratch.size() >= 2)
     {
       _rs.push_back(scratch[0]);
       _xs.push_back(scratch[1]);
       if (scratch.size() >= 3)
         _ys.push_back(scratch[2]);
       else
         _ys.push_back(0.0);
       if (scratch.size() >= 4)
         _zs.push_back(scratch[3]);
       else
         _zs.push_back(0.0);
     }
   }
  
   file.close();
  
   const int num_pts = _zs.size();
   _bottom_point(0) = _xs[num_pts - 1];
   _bottom_point(1) = _ys[num_pts - 1];
   _bottom_point(2) = _zs[num_pts - 1];
  
   // construct the line-segment lengths between each point
   _half_seg_len.resize(std::max(num_pts - 1, 1));
   for (unsigned int i = 0; i + 1 < _xs.size(); ++i)
   {
     _half_seg_len[i] =
         0.5 * std::sqrt(std::pow(_xs[i + 1] - _xs[i], 2) + std::pow(_ys[i + 1] - _ys[i], 2) +
                         std::pow(_zs[i + 1] - _zs[i], 2));
     if (_half_seg_len[i] == 0)
       mooseError("Peaceman-type borehole has a zero-segment length at (x,y,z) = ",
                  _xs[i],
                  " ",
                  _ys[i],
                  " ",
                  _zs[i],
                  "\n");
   }
   if (num_pts == 1)
     _half_seg_len[0] = _borehole_length;
  
   // construct the rotation matrix needed to rotate the permeability
   _rot_matrix.resize(std::max(num_pts - 1, 1));
   for (unsigned int i = 0; i + 1 < _xs.size(); ++i)
   {
     const RealVectorValue v2(_xs[i + 1] - _xs[i], _ys[i + 1] - _ys[i], _zs[i + 1] - _zs[i]);
     _rot_matrix[i] = RotationMatrix::rotVecToZ(v2);
   }
   if (num_pts == 1)
     _rot_matrix[0] = RotationMatrix::rotVecToZ(_borehole_direction);
 }

Member Function Documentation

◆ addPoints()

void PeacemanBorehole::addPoints ( )

protectedvirtual

Add Dirac Points to the borehole.

Definition at line 180 of file PeacemanBorehole.C.

 {
   // This function gets called just before the DiracKernel is evaluated
   // so this is a handy place to zero this out.
   _total_outflow_mass.zero();
  
   // Add point using the unique ID "i", let the DiracKernel take
   // care of the caching.  This should be fast after the first call,
   // as long as the points don't move around.
   for (unsigned int i = 0; i < _zs.size(); i++)
     addPoint(Point(_xs[i], _ys[i], _zs[i]), i);
 }

◆ parseNextLineReals()

bool PeacemanBorehole::parseNextLineReals	(	std::ifstream &	ifs,
		std::vector< Real > &	myvec
	)

protected

reads a space-separated line of floats from ifs and puts in myvec

Definition at line 158 of file PeacemanBorehole.C.

 {
   std::string line;
   myvec.clear();
   bool gotline(false);
   if (getline(ifs, line))
   {
     gotline = true;
  
     // Harvest floats separated by whitespace
     std::istringstream iss(line);
     Real f;
     while (iss >> f)
     {
       myvec.push_back(f);
     }
   }
   return gotline;
 }

Referenced by PeacemanBorehole().

◆ wellConstant()

Real PeacemanBorehole::wellConstant	(	const RealTensorValue &	perm,
		const RealTensorValue &	rot,
		const Real &	half_len,
		const Elem *	ele,
		const Real &	rad
	)

protected

Calculates Peaceman's form of the borehole well constant Z Chen, Y Zhang, Well flow models for various numerical methods, Int J Num Analysis and Modeling, 3 (2008) 375-388.

Definition at line 194 of file PeacemanBorehole.C.

 {
   if (_well_constant > 0)
     return _well_constant;
  
   // rot_perm has its "2" component lying along the half segment
   // we want to determine the eigenvectors of rot(0:1, 0:1), since, when
   // rotated back to the original frame we will determine the element
   // lengths along these directions
   const RealTensorValue rot_perm = (rot * perm) * rot.transpose();
   const Real trace2D = rot_perm(0, 0) + rot_perm(1, 1);
   const Real det2D = rot_perm(0, 0) * rot_perm(1, 1) - rot_perm(0, 1) * rot_perm(1, 0);
   const Real sq = std::sqrt(std::max(0.25 * trace2D * trace2D - det2D,
                                      0.0)); // the std::max accounts for wierdo precision loss
   const Real eig_val1 = 0.5 * trace2D + sq;
   const Real eig_val2 = 0.5 * trace2D - sq;
   RealVectorValue eig_vec1, eig_vec2;
   if (sq > std::abs(trace2D) * 1E-7) // matrix is not a multiple of the identity (1E-7 accounts for
                                      // precision in a crude way)
   {
     if (rot_perm(1, 0) != 0)
     {
       eig_vec1(0) = eig_val1 - rot_perm(1, 1);
       eig_vec1(1) = rot_perm(1, 0);
       eig_vec2(0) = eig_val2 - rot_perm(1, 1);
       eig_vec2(1) = rot_perm(1, 0);
     }
     else if (rot_perm(0, 1) != 0)
     {
       eig_vec1(0) = rot_perm(0, 1);
       eig_vec1(1) = eig_val1 - rot_perm(0, 0);
       eig_vec2(0) = rot_perm(0, 1);
       eig_vec2(1) = eig_val2 - rot_perm(0, 0);
     }
     else // off diagonal terms are both zero
     {
       eig_vec1(0) = 1;
       eig_vec2(1) = 1;
     }
   }
   else // matrix is basically a multiple of the identity
   {
     eig_vec1(0) = 1;
     eig_vec2(1) = 1;
   }
  
   // finally, rotate these to original frame and normalise
   eig_vec1 = rot.transpose() * eig_vec1;
   eig_vec1 /= std::sqrt(eig_vec1 * eig_vec1);
   eig_vec2 = rot.transpose() * eig_vec2;
   eig_vec2 /= std::sqrt(eig_vec2 * eig_vec2);
  
   // find the "length" of the element in these directions
   // TODO - probably better to use variance than max&min
   Real max1 = eig_vec1 * ele->point(0);
   Real max2 = eig_vec2 * ele->point(0);
   Real min1 = max1;
   Real min2 = max2;
   Real proj;
   for (unsigned int i = 1; i < ele->n_nodes(); i++)
   {
     proj = eig_vec1 * ele->point(i);
     max1 = (max1 < proj) ? proj : max1;
     min1 = (min1 < proj) ? min1 : proj;
  
     proj = eig_vec2 * ele->point(i);
     max2 = (max2 < proj) ? proj : max2;
     min2 = (min2 < proj) ? min2 : proj;
   }
   const Real ll1 = max1 - min1;
   const Real ll2 = max2 - min2;
  
   const Real r0 = _re_constant * std::sqrt(std::sqrt(eig_val1 / eig_val2) * std::pow(ll2, 2) +
                                            std::sqrt(eig_val2 / eig_val1) * std::pow(ll1, 2)) /
                   (std::pow(eig_val1 / eig_val2, 0.25) + std::pow(eig_val2 / eig_val1, 0.25));
  
   const Real effective_perm = std::sqrt(det2D);
  
   const Real halfPi = acos(0.0);
  
   if (r0 <= rad)
     mooseError("The effective element size (about 0.2-times-true-ele-size) for an element "
                "containing a Peaceman-type borehole must be (much) larger than the borehole radius "
                "for the Peaceman formulation to be correct.  Your element has effective size ",
                r0,
                " and the borehole radius is ",
                rad,
                "\n");
  
   return 4 * halfPi * effective_perm * half_len / std::log(r0 / rad);
 }

Referenced by Q2PBorehole::computeQpResidual(), RichardsBorehole::computeQpResidual(), Q2PBorehole::jac(), and RichardsBorehole::jac().

Member Data Documentation

◆ _borehole_direction

const RealVectorValue PeacemanBorehole::_borehole_direction

private

borehole direction. Note this is only used if there is only one borehole point

Definition at line 49 of file PeacemanBorehole.h.

Referenced by PeacemanBorehole().

◆ _borehole_length

const Real PeacemanBorehole::_borehole_length

private

borehole length. Note this is only used if there is only one borehole point

Definition at line 46 of file PeacemanBorehole.h.

Referenced by PeacemanBorehole().

◆ _bottom_point

Point PeacemanBorehole::_bottom_point

protected

the bottom point of the borehole (where bottom_pressure is defined)

Definition at line 94 of file PeacemanBorehole.h.

Referenced by Q2PBorehole::computeQpResidual(), RichardsBorehole::computeQpResidual(), Q2PBorehole::jac(), RichardsBorehole::jac(), and PeacemanBorehole().

◆ _character

const Function& PeacemanBorehole::_character

protected

If positive then the borehole acts as a sink (producion well) for porepressure > borehole pressure, and does nothing otherwise If negative then the borehole acts as a source (injection well) for porepressure < borehole pressure, and does nothing otherwise The flow rate to/from the borehole is multiplied by |character|, so usually character = +/- 1.

Definition at line 66 of file PeacemanBorehole.h.

Referenced by RichardsBorehole::computeQpJacobian(), Q2PBorehole::computeQpResidual(), RichardsBorehole::computeQpResidual(), Q2PBorehole::jac(), and RichardsBorehole::jac().

◆ _half_seg_len

std::vector<Real> PeacemanBorehole::_half_seg_len

protected

0.5*(length of polyline segments between points)

Definition at line 97 of file PeacemanBorehole.h.

Referenced by Q2PBorehole::computeQpResidual(), RichardsBorehole::computeQpResidual(), Q2PBorehole::jac(), RichardsBorehole::jac(), and PeacemanBorehole().

◆ _p_bot

const Real PeacemanBorehole::_p_bot

protected

bottomhole pressure of borehole

Definition at line 69 of file PeacemanBorehole.h.

Referenced by Q2PBorehole::computeQpResidual(), RichardsBorehole::computeQpResidual(), Q2PBorehole::jac(), and RichardsBorehole::jac().

◆ _point_file

const std::string PeacemanBorehole::_point_file

private

File defining the geometry of the borehole.

Each row has format radius x y z and the list of such points defines a polyline that is the borehole

Definition at line 56 of file PeacemanBorehole.h.

Referenced by PeacemanBorehole().

◆ _re_constant

const Real PeacemanBorehole::_re_constant

private

borehole constant

Definition at line 40 of file PeacemanBorehole.h.

Referenced by wellConstant().

◆ _rot_matrix

std::vector<RealTensorValue> PeacemanBorehole::_rot_matrix

protected

rotation matrix used in well_constant calculation

Definition at line 100 of file PeacemanBorehole.h.

Referenced by Q2PBorehole::computeQpResidual(), RichardsBorehole::computeQpResidual(), Q2PBorehole::jac(), RichardsBorehole::jac(), and PeacemanBorehole().

◆ _rs

std::vector<Real> PeacemanBorehole::_rs

protected

radii of the borehole

Definition at line 82 of file PeacemanBorehole.h.

Referenced by Q2PBorehole::computeQpResidual(), RichardsBorehole::computeQpResidual(), Q2PBorehole::jac(), RichardsBorehole::jac(), and PeacemanBorehole().

◆ _total_outflow_mass

RichardsSumQuantity& PeacemanBorehole::_total_outflow_mass

protected

This is used to hold the total fluid flowing into the borehole Hence, it is positive for production wells where fluid is flowing from porespace into the borehole and removed from the model.

Definition at line 79 of file PeacemanBorehole.h.

Referenced by addPoints(), Q2PBorehole::computeQpResidual(), RichardsBorehole::computeQpResidual(), and PeacemanBorehole().

◆ _unit_weight

const RealVectorValue PeacemanBorehole::_unit_weight

protected

unit weight of fluid in borehole (for calculating bottomhole pressure at each Dirac Point)

Definition at line 72 of file PeacemanBorehole.h.

Referenced by Q2PBorehole::computeQpResidual(), RichardsBorehole::computeQpResidual(), Q2PBorehole::jac(), and RichardsBorehole::jac().

◆ _well_constant

const Real PeacemanBorehole::_well_constant

private

well constant

Definition at line 43 of file PeacemanBorehole.h.

Referenced by wellConstant().

◆ _xs

std::vector<Real> PeacemanBorehole::_xs

protected

x points of the borehole

Definition at line 85 of file PeacemanBorehole.h.

Referenced by addPoints(), and PeacemanBorehole().

◆ _ys

std::vector<Real> PeacemanBorehole::_ys

protected

y points of the borehole

Definition at line 88 of file PeacemanBorehole.h.

Referenced by addPoints(), and PeacemanBorehole().

◆ _zs

std::vector<Real> PeacemanBorehole::_zs

protected

z points of borehole

Definition at line 91 of file PeacemanBorehole.h.

Referenced by addPoints(), Q2PBorehole::computeQpResidual(), RichardsBorehole::computeQpResidual(), Q2PBorehole::jac(), RichardsBorehole::jac(), and PeacemanBorehole().

The documentation for this class was generated from the following files:

richards/include/dirackernels/PeacemanBorehole.h
richards/src/dirackernels/PeacemanBorehole.C

Public Member Functions

Protected Member Functions

Protected Attributes

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ PeacemanBorehole()

Member Function Documentation

◆ addPoints()

◆ parseNextLineReals()

◆ wellConstant()

Member Data Documentation

◆ _borehole_direction

◆ _borehole_length

◆ _bottom_point

◆ _character

◆ _half_seg_len

◆ _p_bot

◆ _point_file

◆ _re_constant

◆ _rot_matrix

◆ _rs

◆ _total_outflow_mass

◆ _unit_weight

◆ _well_constant

◆ _xs

◆ _ys

◆ _zs