mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Added hydrostatic correction for vfp table depth
Closer to reproducing proper results, but som oscillating behaviour in plots of actual BHP.
This commit is contained in:
babrodtk
parent
748440eea6
commit
657a7c58b8
@ -260,6 +260,7 @@ namespace Opm {
|
||||
V isRs_;
|
||||
V isRv_;
|
||||
V isSg_;
|
||||
V well_perforation_densities_; //Density of each well perforation
|
||||
V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
|
||||
|
||||
LinearisedBlackoilResidual residual_;
|
||||
|
||||
@ -657,7 +657,19 @@ namespace detail {
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
namespace detail {
|
||||
double getGravity(const double* g, const int dim) {
|
||||
double grav = 0.0;
|
||||
if (g) {
|
||||
// Guard against gravity in anything but last dimension.
|
||||
for (int dd = 0; dd < dim - 1; ++dd) {
|
||||
assert(g[dd] == 0.0);
|
||||
}
|
||||
grav = g[dim - 1];
|
||||
}
|
||||
return grav;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -731,22 +743,21 @@ namespace detail {
|
||||
// Surface density.
|
||||
std::vector<double> surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases);
|
||||
// Gravity
|
||||
double grav = 0.0;
|
||||
const double* g = geo_.gravity();
|
||||
const int dim = dimensions(grid_);
|
||||
if (g) {
|
||||
// Guard against gravity in anything but last dimension.
|
||||
for (int dd = 0; dd < dim - 1; ++dd) {
|
||||
assert(g[dd] == 0.0);
|
||||
}
|
||||
grav = g[dim - 1];
|
||||
}
|
||||
double grav = detail::getGravity(geo_.gravity(), dimensions(grid_));
|
||||
|
||||
// 2. Compute pressure deltas, and store the results.
|
||||
std::vector<double> cdp = WellDensitySegmented
|
||||
::computeConnectionPressureDelta(wells(), xw, fluid_.phaseUsage(),
|
||||
b_perf, rsmax_perf, rvmax_perf, perf_depth,
|
||||
surf_dens, grav);
|
||||
// 2. Compute densities
|
||||
std::vector<double> cd =
|
||||
WellDensitySegmented::computeConnectionDensities(
|
||||
wells(), xw, fluid_.phaseUsage(),
|
||||
b_perf, rsmax_perf, rvmax_perf, surf_dens);
|
||||
|
||||
// 3. Compute pressure deltas
|
||||
std::vector<double> cdp =
|
||||
WellDensitySegmented::computeConnectionPressureDelta(
|
||||
wells(), perf_depth, cd, grav);
|
||||
|
||||
// 4. Store the results
|
||||
well_perforation_densities_ = Eigen::Map<const V>(cd.data(), nperf);
|
||||
well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf);
|
||||
}
|
||||
|
||||
@ -1166,6 +1177,25 @@ namespace detail {
|
||||
} // namespace detail
|
||||
|
||||
|
||||
namespace detail {
|
||||
double computeHydrostaticCorrection(const Wells& wells, const int w, const double vfp_ref_depth,
|
||||
const ADB::V& well_perforation_densities, const double gravity) {
|
||||
//For the initial iteration, we have no perforation densities.
|
||||
if (well_perforation_densities.size() > w) {
|
||||
const double well_ref_depth = wells.depth_ref[w];
|
||||
const double dh = vfp_ref_depth - well_ref_depth;
|
||||
const int perf = wells.well_connpos[w];
|
||||
const double rho = well_perforation_densities[perf];
|
||||
const double dp = rho*gravity*dh;
|
||||
|
||||
return dp;
|
||||
}
|
||||
else {
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
} //Namespace
|
||||
|
||||
|
||||
|
||||
@ -1218,6 +1248,9 @@ namespace detail {
|
||||
current = xw.currentControls()[w];
|
||||
}
|
||||
|
||||
//Get gravity for THP hydrostatic corrrection
|
||||
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
|
||||
|
||||
// Updating well state and primary variables.
|
||||
// Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
|
||||
const double target = well_controls_iget_target(wc, current);
|
||||
@ -1248,11 +1281,24 @@ namespace detail {
|
||||
|
||||
//Set *BHP* target by calculating bhp from THP
|
||||
const WellType& well_type = wells().type[w];
|
||||
if (well_type == PRODUCER) {
|
||||
xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq);
|
||||
|
||||
//Gather variables for hydrostatic reconstruction
|
||||
|
||||
if (well_type == INJECTOR) {
|
||||
double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth();
|
||||
double dp = detail::computeHydrostaticCorrection(
|
||||
wells(), w, vfp_ref_depth,
|
||||
well_perforation_densities_, gravity);
|
||||
|
||||
xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
|
||||
}
|
||||
else if (well_type == INJECTOR) {
|
||||
xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp);
|
||||
else if (well_type == PRODUCER) {
|
||||
double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth();
|
||||
double dp = detail::computeHydrostaticCorrection(
|
||||
wells(), w, vfp_ref_depth,
|
||||
well_perforation_densities_, gravity);
|
||||
|
||||
xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
|
||||
}
|
||||
else {
|
||||
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
|
||||
@ -1386,10 +1432,14 @@ namespace detail {
|
||||
//THP calculation variables
|
||||
std::vector<int> inj_table_id(nw, -1);
|
||||
std::vector<int> prod_table_id(nw, -1);
|
||||
ADB::V thp_inj_v = ADB::V::Zero(nw);
|
||||
ADB::V thp_prod_v = ADB::V::Zero(nw);
|
||||
ADB::V thp_inj_target_v = ADB::V::Zero(nw);
|
||||
ADB::V thp_prod_target_v = ADB::V::Zero(nw);
|
||||
ADB::V alq_v = ADB::V::Zero(nw);
|
||||
|
||||
//Hydrostatic correction variables
|
||||
ADB::V rho_v = ADB::V::Zero(nw);
|
||||
ADB::V vfp_ref_depth_v = ADB::V::Zero(nw);
|
||||
|
||||
//Target vars
|
||||
ADB::V bhp_targets = ADB::V::Zero(nw);
|
||||
ADB::V rate_targets = ADB::V::Zero(nw);
|
||||
@ -1422,19 +1472,29 @@ namespace detail {
|
||||
|
||||
case THP:
|
||||
{
|
||||
const int perf = wells().well_connpos[w];
|
||||
rho_v[w] = well_perforation_densities_[perf];
|
||||
|
||||
const int table_id = well_controls_iget_vfp(wc, current);
|
||||
const double target = well_controls_iget_target(wc, current);
|
||||
|
||||
const WellType& well_type = wells().type[w];
|
||||
if (well_type == INJECTOR) {
|
||||
inj_table_id[w] = well_controls_iget_vfp(wc, current);
|
||||
thp_inj_v[w] = well_controls_iget_target(wc, current);
|
||||
inj_table_id[w] = table_id;
|
||||
thp_inj_target_v[w] = target;
|
||||
alq_v[w] = -1e100;
|
||||
|
||||
vfp_ref_depth_v[w] = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
|
||||
|
||||
thp_inj_elems.push_back(w);
|
||||
}
|
||||
else if (well_type == PRODUCER) {
|
||||
prod_table_id[w] = well_controls_iget_vfp(wc, current);
|
||||
thp_prod_v[w] = well_controls_iget_target(wc, current);
|
||||
prod_table_id[w] = table_id;
|
||||
thp_prod_target_v[w] = target;
|
||||
alq_v[w] = well_controls_iget_alq(wc, current);
|
||||
|
||||
vfp_ref_depth_v[w] = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
|
||||
|
||||
thp_prod_elems.push_back(w);
|
||||
}
|
||||
else {
|
||||
@ -1468,23 +1528,33 @@ namespace detail {
|
||||
}
|
||||
|
||||
//Calculate BHP target from THP
|
||||
const ADB thp_inj = ADB::constant(thp_inj_v);
|
||||
const ADB thp_prod = ADB::constant(thp_prod_v);
|
||||
const ADB thp_inj_target = ADB::constant(thp_inj_target_v);
|
||||
const ADB thp_prod_target = ADB::constant(thp_prod_target_v);
|
||||
const ADB alq = ADB::constant(alq_v);
|
||||
const ADB thp_inj_targets = vfp_properties_->getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj);
|
||||
const ADB thp_prod_targets = vfp_properties_->getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod, alq);
|
||||
const ADB bhp_from_thp_inj = vfp_properties_->getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj_target);
|
||||
const ADB bhp_from_thp_prod = vfp_properties_->getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq);
|
||||
|
||||
//Perform hydrostatic correction to computed targets
|
||||
const ADB well_ref_depth = ADB::constant(ADB::V::Map(wells().depth_ref, nw));
|
||||
const ADB vfp_ref_depth = ADB::constant(vfp_ref_depth_v);
|
||||
const ADB dh = vfp_ref_depth - well_ref_depth;
|
||||
const ADB rho = ADB::constant(rho_v);
|
||||
double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
|
||||
const ADB dp = rho*gravity*dh;
|
||||
const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw);
|
||||
const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw);
|
||||
|
||||
//Calculate residuals
|
||||
const ADB thp_inj_residual = state.bhp - thp_inj_targets;
|
||||
const ADB thp_prod_residual = state.bhp - thp_prod_targets;
|
||||
const ADB thp_inj_residual = state.bhp - bhp_from_thp_inj + dp_inj;
|
||||
const ADB thp_prod_residual = state.bhp - bhp_from_thp_prod + dp_prod;
|
||||
const ADB bhp_residual = state.bhp - bhp_targets;
|
||||
const ADB rate_residual = rate_distr * state.qs - rate_targets;
|
||||
|
||||
//Select the right residual for each well
|
||||
residual_.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, bhp_residual.size()) +
|
||||
superset(subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, thp_inj_residual.size()) +
|
||||
superset(subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, thp_prod_residual.size()) +
|
||||
superset(subset(rate_residual, rate_elems), rate_elems, rate_residual.size());
|
||||
residual_.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, nw) +
|
||||
superset(subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, nw) +
|
||||
superset(subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, nw) +
|
||||
superset(subset(rate_residual, rate_elems), rate_elems, nw);
|
||||
|
||||
// For wells that are dead (not flowing), and therefore not communicating
|
||||
// with the reservoir, we set the equation to be equal to the well's total
|
||||
|
||||
@ -94,7 +94,6 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
|
||||
|
||||
|
||||
|
||||
|
||||
VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
|
||||
const ADB& aqua,
|
||||
const ADB& liquid,
|
||||
@ -118,7 +117,7 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
|
||||
//Get the table for each well
|
||||
std::vector<const VFPInjTable*> well_tables(nw, NULL);
|
||||
for (int i=0; i<nw; ++i) {
|
||||
if (table_id[i] >= 0) {
|
||||
if (table_id[i] > 0) {
|
||||
well_tables[i] = detail::getTable(m_tables, table_id[i]);
|
||||
}
|
||||
}
|
||||
@ -226,4 +225,14 @@ double VFPInjProperties::thp(int table_id,
|
||||
|
||||
|
||||
|
||||
const VFPInjTable* VFPInjProperties::getTable(const int table_id) const {
|
||||
return detail::getTable(m_tables, table_id);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
} //Namespace Opm
|
||||
|
||||
@ -126,6 +126,12 @@ public:
|
||||
const double& vapour,
|
||||
const double& bhp) const;
|
||||
|
||||
/**
|
||||
* Returns the table associated with the ID, or throws an exception if
|
||||
* the table does not exist
|
||||
*/
|
||||
const VFPInjTable* getTable(const int table_id) const;
|
||||
|
||||
private:
|
||||
// Map which connects the table number with the table itself
|
||||
std::map<int, const VFPInjTable*> m_tables;
|
||||
|
||||
@ -235,6 +235,14 @@ double VFPProdProperties::thp(int table_id,
|
||||
|
||||
|
||||
|
||||
const VFPProdTable* VFPProdProperties::getTable(const int table_id) const {
|
||||
return detail::getTable(m_tables, table_id);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
@ -137,6 +137,11 @@ public:
|
||||
const double& bhp,
|
||||
const double& alq) const;
|
||||
|
||||
/**
|
||||
* Returns the table associated with the ID, or throws an exception if
|
||||
* the table does not exist
|
||||
*/
|
||||
const VFPProdTable* getTable(const int table_id) const;
|
||||
|
||||
private:
|
||||
// Map which connects the table number with the table itself
|
||||
|
||||
@ -26,16 +26,17 @@
|
||||
#include <cmath>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
std::vector<double>
|
||||
Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
|
||||
Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
|
||||
const WellStateFullyImplicitBlackoil& wstate,
|
||||
const PhaseUsage& phase_usage,
|
||||
const std::vector<double>& b_perf,
|
||||
const std::vector<double>& rsmax_perf,
|
||||
const std::vector<double>& rvmax_perf,
|
||||
const std::vector<double>& z_perf,
|
||||
const std::vector<double>& surf_dens,
|
||||
const double gravity)
|
||||
const std::vector<double>& surf_dens)
|
||||
{
|
||||
// Verify that we have consistent input.
|
||||
const int np = wells.number_of_phases;
|
||||
@ -53,9 +54,6 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
|
||||
if (nperf*np != int(b_perf.size())) {
|
||||
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. b_perf.");
|
||||
}
|
||||
if (nperf != int(z_perf.size())) {
|
||||
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. z_perf.");
|
||||
}
|
||||
if ((!rsmax_perf.empty()) || (!rvmax_perf.empty())) {
|
||||
// Need both oil and gas phases.
|
||||
if (!phase_usage.phase_used[BlackoilPhases::Liquid]) {
|
||||
@ -66,14 +64,6 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
|
||||
}
|
||||
}
|
||||
|
||||
// Algorithm:
|
||||
|
||||
// We'll assume the perforations are given in order from top to
|
||||
// bottom for each well. By top and bottom we do not necessarily
|
||||
// mean in a geometric sense (depth), but in a topological sense:
|
||||
// the 'top' perforation is nearest to the surface topologically.
|
||||
// Our goal is to compute a pressure delta for each perforation.
|
||||
|
||||
// 1. Compute the flow (in surface volume units for each
|
||||
// component) exiting up the wellbore from each perforation,
|
||||
// taking into account flow from lower in the well, and
|
||||
@ -145,7 +135,38 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
|
||||
}
|
||||
}
|
||||
|
||||
// 3. Compute pressure differences between perforations.
|
||||
return dens;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
std::vector<double>
|
||||
Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
|
||||
const std::vector<double>& z_perf,
|
||||
const std::vector<double>& dens_perf,
|
||||
const double gravity) {
|
||||
const int np = wells.number_of_phases;
|
||||
const int nw = wells.number_of_wells;
|
||||
const int nperf = wells.well_connpos[nw];
|
||||
|
||||
if (nperf != int(z_perf.size())) {
|
||||
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. z_perf.");
|
||||
}
|
||||
if (nperf != int(dens_perf.size())) {
|
||||
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. dens_perf.");
|
||||
}
|
||||
|
||||
// Algorithm:
|
||||
|
||||
// We'll assume the perforations are given in order from top to
|
||||
// bottom for each well. By top and bottom we do not necessarily
|
||||
// mean in a geometric sense (depth), but in a topological sense:
|
||||
// the 'top' perforation is nearest to the surface topologically.
|
||||
// Our goal is to compute a pressure delta for each perforation.
|
||||
|
||||
// 1. Compute pressure differences between perforations.
|
||||
// dp_perf will contain the pressure difference between a
|
||||
// perforation and the one above it, except for the first
|
||||
// perforation for each well, for which it will be the
|
||||
@ -155,11 +176,12 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
|
||||
for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w+1]; ++perf) {
|
||||
const double z_above = perf == wells.well_connpos[w] ? wells.depth_ref[w] : z_perf[perf - 1];
|
||||
const double dz = z_perf[perf] - z_above;
|
||||
dp_perf[perf] = dz * dens[perf] * gravity;
|
||||
dp_perf[perf] = dz * dens_perf[perf] * gravity;
|
||||
//dens[wells.well_connpos[w]]
|
||||
}
|
||||
}
|
||||
|
||||
// 4. Compute pressure differences to the reference point (bhp) by
|
||||
// 2. Compute pressure differences to the reference point (bhp) by
|
||||
// accumulating the already computed adjacent pressure
|
||||
// differences, storing the result in dp_perf.
|
||||
// This accumulation must be done per well.
|
||||
|
||||
@ -39,7 +39,7 @@ namespace Opm
|
||||
class WellDensitySegmented
|
||||
{
|
||||
public:
|
||||
/// Compute pressure deltas.
|
||||
/// Compute well segment densities
|
||||
/// Notation: N = number of perforations, P = number of phases.
|
||||
/// \param[in] wells struct with static well info
|
||||
/// \param[in] wstate dynamic well solution information, only perfRates() is used
|
||||
@ -47,17 +47,24 @@ namespace Opm
|
||||
/// \param[in] b_perf inverse ('little b') formation volume factor, size NP, P values per perforation
|
||||
/// \param[in] rsmax_perf saturation point for rs (gas in oil) at each perforation, size N
|
||||
/// \param[in] rvmax_perf saturation point for rv (oil in gas) at each perforation, size N
|
||||
/// \param[in] z_perf depth values for each perforation, size N
|
||||
/// \param[in] surf_dens surface densities for active components, size P
|
||||
/// \param[in] gravity gravity acceleration constant
|
||||
static std::vector<double> computeConnectionPressureDelta(const Wells& wells,
|
||||
static std::vector<double> computeConnectionDensities(const Wells& wells,
|
||||
const WellStateFullyImplicitBlackoil& wstate,
|
||||
const PhaseUsage& phase_usage,
|
||||
const std::vector<double>& b_perf,
|
||||
const std::vector<double>& rsmax_perf,
|
||||
const std::vector<double>& rvmax_perf,
|
||||
const std::vector<double>& surf_dens);
|
||||
|
||||
/// Compute pressure deltas.
|
||||
/// Notation: N = number of perforations, P = number of phases.
|
||||
/// \param[in] wells struct with static well info
|
||||
/// \param[in] z_perf depth values for each perforation, size N
|
||||
/// \param[in] dens_perf densities for each perforation, size N (typically computed using computeConnectionDensities)
|
||||
/// \param[in] gravity gravity acceleration constant
|
||||
static std::vector<double> computeConnectionPressureDelta(const Wells& wells,
|
||||
const std::vector<double>& z_perf,
|
||||
const std::vector<double>& surf_dens,
|
||||
const std::vector<double>& dens_perf,
|
||||
const double gravity);
|
||||
};
|
||||
|
||||
|
||||
@ -90,7 +90,15 @@ BOOST_AUTO_TEST_CASE(TestPressureDeltas)
|
||||
const std::vector<double> z_perf = { 10, 30, 50, 70, 90, 10, 30, 50, 70, 90 };
|
||||
const std::vector<double> surf_dens = { 1000.0, 800.0, 10.0 };
|
||||
const double gravity = Opm::unit::gravity;
|
||||
const std::vector<double> dp = WellDensitySegmented::computeConnectionPressureDelta(*wells, wellstate, pu, b_perf, rsmax_perf, rvmax_perf, z_perf, surf_dens, gravity);
|
||||
|
||||
std::vector<double> cd =
|
||||
WellDensitySegmented::computeConnectionDensities(
|
||||
*wells, wellstate, pu,
|
||||
b_perf, rsmax_perf, rvmax_perf, surf_dens);
|
||||
std::vector<double> dp =
|
||||
WellDensitySegmented::computeConnectionPressureDelta(
|
||||
*wells, z_perf, cd, gravity);
|
||||
|
||||
const std::vector<double> answer = { 20e3*gravity, 62e3*gravity, 106e3*gravity, 152e3*gravity, 200e3*gravity,
|
||||
20e3*gravity, 62e3*gravity, 106e3*gravity, 152e3*gravity, 200e3*gravity };
|
||||
BOOST_REQUIRE_EQUAL(dp.size(), answer.size());
|
||||
|
||||
Loading…
Reference in New Issue
Block a user