Use Well CF and Kh calculations from opm-common

This commit is contained in:
Joakim Hove 2018-08-31 08:50:44 +02:00
parent 25248066d8
commit b9df94a70f
5 changed files with 14 additions and 180 deletions

View File

@ -59,7 +59,7 @@ namespace Opm
const WellConnections& completion_set = well_ecl_->getConnections(current_step_);
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const Connection& connection = completion_set.get(perf);
const int segment_index = segmentNumberToIndex(connection.segment_number);
const int segment_index = segmentNumberToIndex(connection.segment());
segment_perforations_[segment_index].push_back(perf);
}
@ -80,7 +80,7 @@ namespace Opm
for (int seg = 0; seg < numberOfSegments(); ++seg) {
const double segment_depth = segmentSet()[seg].depth();
for (const int perf : segment_perforations_[seg]) {
perf_depth_[perf] = completion_set.get(perf).center_depth;
perf_depth_[perf] = completion_set.get(perf).depth();
perforation_segment_depth_diffs_[perf] = perf_depth_[perf] - segment_depth;
}
}

View File

@ -552,7 +552,7 @@ namespace Opm
for (const auto& completion : completions) {
int complnum = completion.first;
for (int perf = 0; perf < perf_number; ++perf) {
if (complnum == connections.get ( perf ).complnum) {
if (complnum == connections.get ( perf ).complnum()) {
water_cut_in_completions[complnumIdx] += water_cut_perf[perf];
}
}
@ -730,7 +730,7 @@ namespace Opm
bool allCompletionsClosed = true;
const auto& connections = well_ecl_->getConnections(current_step_);
for (const auto& connection : connections) {
if (!wellTestState.hasCompletion(name(), connection.complnum)) {
if (!wellTestState.hasCompletion(name(), connection.complnum())) {
allCompletionsClosed = false;
}
}
@ -802,7 +802,7 @@ namespace Opm
const auto& connectionSet = well_ecl_->getConnections(current_step_);
for (size_t c=0; c<connectionSet.size(); c++) {
const auto& connection = connectionSet.get(c);
if (connection.state == WellCompletion::OPEN) {
if (connection.state() == WellCompletion::OPEN) {
const int i = connection.getI();
const int j = connection.getJ();
const int k = connection.getK();
@ -817,19 +817,14 @@ namespace Opm
const int cell = cgit->second;
{
double radius = 0.5*connection.getDiameter();
if (radius <= 0.0) {
radius = 0.5*unit::feet;
OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
}
double radius = connection.rw();
const std::array<double, 3> cubical =
WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
double re; // area equivalent radius of the grid block
double perf_length; // the length of the well perforation
switch (connection.dir) {
switch (connection.dir()) {
case Opm::WellCompletion::DirectionEnum::X:
re = std::sqrt(cubical[1] * cubical[2] / M_PI);
perf_length = cubical[0];
@ -918,7 +913,7 @@ namespace Opm
const auto& connections = well_ecl_->getConnections(current_step_);
int perfIdx = 0;
for (const auto& connection : connections) {
if (wellTestState.hasCompletion(name(), connection.complnum)) {
if (wellTestState.hasCompletion(name(), connection.complnum())) {
well_index_[perfIdx] = 0.0;
}
perfIdx++;

View File

@ -559,7 +559,7 @@ namespace Opm
std::vector<std::vector<int>> segment_perforations(well_nseg);
for (int perf = 0; perf < nperf; ++perf) {
const Connection& connection = completion_set.get(perf);
const int segment_index = segment_set.segmentNumberToIndex(connection.segment_number);
const int segment_index = segment_set.segmentNumberToIndex(connection.segment());
segment_perforations[segment_index].push_back(perf);
}

View File

@ -175,134 +175,6 @@ namespace WellsManagerDetail
} // namespace InjectionControl
// Compute direction permutation corresponding to completion's
// direction. First two elements of return value are directions
// perpendicular to completion while last element is direction
// along completion.
inline std::array< std::array<double,3>::size_type, 3 >
directionIndices(const Opm::WellCompletion::DirectionEnum direction)
{
typedef std::array<double,3>::size_type idx_t;
typedef std::array<idx_t,3> permutation;
switch (direction) {
case Opm::WellCompletion::DirectionEnum::X:
return permutation {{ idx_t(1), idx_t(2), idx_t(0) }};
case Opm::WellCompletion::DirectionEnum::Y:
return permutation {{ idx_t(2), idx_t(0), idx_t(1) }};
case Opm::WellCompletion::DirectionEnum::Z:
return permutation {{ idx_t(0), idx_t(1), idx_t(2) }};
}
// All enum values should be handled above. Therefore
// we should never reach this one. Anyway for the sake
// of reduced warnings we throw an exception.
throw std::invalid_argument("unhandled enum value");
}
// Permute (diagonal) permeability components according to
// completion's direction.
inline std::array<double,3>
permComponents(const Opm::WellCompletion::DirectionEnum direction,
const double* perm)
{
const auto p = directionIndices(direction);
const std::array<double,3>::size_type d = 3;
std::array<double,3>
K = {{ perm[ p[0]*(d + 1) ] ,
perm[ p[1]*(d + 1) ] ,
perm[ p[2]*(d + 1) ] }};
return K;
}
// Permute cell's geometric extent according to completion's
// direction. Honour net-to-gross ratio.
//
// Note: 'extent' is intentionally accepted by modifiable value
// rather than reference-to-const to support NTG manipulation.
inline std::array<double,3>
effectiveExtent(const Opm::WellCompletion::DirectionEnum direction,
const double ntg,
std::array<double,3> extent)
{
// Vertical extent affected by net-to-gross ratio.
extent[2] *= ntg;
const auto p = directionIndices(direction);
std::array<double,3>
D = {{ extent[ p[0] ] ,
extent[ p[1] ] ,
extent[ p[2] ] }};
return D;
}
// Compute Peaceman's effective radius of single completion.
inline double
effectiveRadius(const std::array<double,3>& K,
const std::array<double,3>& D)
{
const double K01 = K[0] / K[1];
const double K10 = K[1] / K[0];
const double D0_sq = D[0] * D[0];
const double D1_sq = D[1] * D[1];
const double num = std::sqrt((std::sqrt(K10) * D0_sq) +
(std::sqrt(K01) * D1_sq));
const double den = std::pow(K01, 0.25) + std::pow(K10, 0.25);
// Note: Analytic constant 0.28 derived for infintely sized
// formation with repeating well placement.
return 0.28 * (num / den);
}
// Use the Peaceman well model to compute well indices.
// radius is the radius of the well.
// cubical contains [dx, dy, dz] of the cell.
// (Note that the well model asumes that each cell is a cuboid).
// cell_permeability is the permeability tensor of the given cell.
// returns the well index of the cell.
double
computeWellIndex(const double radius,
const std::array<double, 3>& cubical,
const double* cell_permeability,
const double skin_factor,
const Opm::WellCompletion::DirectionEnum direction,
const double ntg)
{
const std::array<double,3>& K =
permComponents(direction, cell_permeability);
const std::array<double,3>& D =
effectiveExtent(direction, ntg, cubical);
const double r0 = effectiveRadius(K, D);
const double Kh = std::sqrt(K[0] * K[1]) * D[2];
// Angle of completion exposed to flow. We assume centre
// placement so there's complete exposure (= 2\pi).
const double angle =
6.2831853071795864769252867665590057683943387987502116419498;
double rw = radius;
if (r0 < rw) {
std::cerr << "Completion radius exceeds effective radius\n"
<< "Reset to effective\n";
rw = r0;
}
// NOTE: The formula is originally derived and valid for
// Cartesian grids only.
return (angle * Kh) / (std::log(r0 / rw) + skin_factor);
}
} // anonymous namespace

View File

@ -54,12 +54,6 @@ namespace WellsManagerDetail
} // namespace InjectionControl
double computeWellIndex(const double radius,
const std::array<double, 3>& cubical,
const double* cell_permeability,
const double skin_factor,
const Opm::WellCompletion::DirectionEnum direction,
const double ntg);
template <int dim, class C2F, class FC>
std::array<double, dim>
@ -165,7 +159,7 @@ void WellsManager::createWellsFromSpecs(std::vector<const Well*>& wells, size_t
// shut completions and open ones stored in this process will have 1 others 0.
for(const auto& completion : well->getConnections(timeStep)) {
if (completion.state == WellCompletion::OPEN) {
if (completion.state() == WellCompletion::OPEN) {
int i = completion.getI();
int j = completion.getJ();
int k = completion.getK();
@ -193,41 +187,14 @@ void WellsManager::createWellsFromSpecs(std::vector<const Well*>& wells, size_t
PerfData pd;
pd.cell = cell;
{
const Value<double>& transmissibilityFactor = completion.getConnectionTransmissibilityFactorAsValueObject();
const double wellPi = completion.wellPi;
if (transmissibilityFactor.hasValue()) {
pd.well_index = transmissibilityFactor.getValue();
} else {
double radius = 0.5*completion.getDiameter();
if (radius <= 0.0) {
radius = 0.5*unit::feet;
OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
}
pd.well_index = completion.CF() * completion.wellPi();
pd.satnumid = completion.satTableId();
std::array<double, 3> cubical =
WellsManagerDetail::getCubeDim<3>(c2f, begin_face_centroids, cell);
// overwrite dz values calculated in getCubeDim.
if (dz.size() > 0) {
cubical[2] = dz[cell];
}
const double* cell_perm = &permeability[dimensions*dimensions*cell];
pd.well_index =
WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm,
completion.getSkinFactor(),
completion.dir,
ntg[cell]);
}
pd.satnumid = completion.sat_tableId;
pd.well_index *= wellPi;
}
wellperf_data[active_well_index].push_back(pd);
}
} else {
if (completion.state != WellCompletion::SHUT) {
OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion.state ) << " not handled");
if (completion.state() != WellCompletion::SHUT) {
OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion.state() ) << " not handled");
}
}
}