From 22f1df612bfbc00263e5677fc4ae5afebc853029 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 3 Oct 2016 09:54:28 +0200 Subject: [PATCH] Allow crossflow if all perforations is crossflowing --- opm/autodiff/StandardWellsDense.hpp | 62 +++++++++++++++++++---------- 1 file changed, 40 insertions(+), 22 deletions(-) diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index f71eab732..5f7409f34 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -226,14 +226,13 @@ namespace Opm { const double volume = 0.002831684659200; // 0.1 cu ft; for (int w = 0; w < nw; ++w) { - double total_flux = 0.0; - + bool allow_cf = allow_cross_flow(w, ebosSimulator); for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { const int cell_idx = wells().well_cells[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); std::vector cq_s(np,0.0); - computeWellFlux(w, wells().WI[perf], intQuants, wellPerforationPressureDiffs()[perf], cq_s); + computeWellFlux(w, wells().WI[perf], intQuants, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); for (int p1 = 0; p1 < np; ++p1) { @@ -244,7 +243,6 @@ namespace Opm { // subtract sum of phase fluxes in the well equations. resWell_[w][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].value; - total_flux += cq_s[p1].value; // assemble the jacobians for (int p2 = 0; p2 < np; ++p2) { @@ -262,23 +260,15 @@ namespace Opm { // Store the perforation pressure for later usage. well_state.perfPress()[perf] = well_state.bhp()[w] + wellPerforationPressureDiffs()[perf]; } - if (std::abs(total_flux) == 0) { // add a trivial equation for the case when there is no flow - for (int p1 = 0; p1 < np; ++p1) { - invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p1)] += 1; - resWell_[w][flowPhaseToEbosCompIdx(p1)] += 0; - } - } - else { - // add vol * dF/dt + Q to the well equations; - for (int p1 = 0; p1 < np; ++p1) { - EvalWell resWell_loc = (wellVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt; - resWell_loc += getQs(w, p1); - for (int p2 = 0; p2 < np; ++p2) { - invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivatives[p2+3]; - } - resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value; + // add vol * dF/dt + Q to the well equations; + for (int p1 = 0; p1 < np; ++p1) { + EvalWell resWell_loc = (wellVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt; + resWell_loc += getQs(w, p1); + for (int p2 = 0; p2 < np; ++p2) { + invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivatives[p2+3]; } + resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value; } } @@ -287,6 +277,34 @@ namespace Opm { } + template + bool allow_cross_flow(const int w, Simulator& ebosSimulator) const { + + if (wells().allow_cf[w]) { + return true; + } + // check for special case where all perforations have cross flow + // then the wells must allow for cross flow + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + const int cell_idx = wells().well_cells[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); + EvalWell bhp = getBhp(w); + // Pressure drawdown (also used to determine direction of flow) + EvalWell well_pressure = bhp + wellPerforationPressureDiffs()[perf]; + EvalWell drawdown = pressure - well_pressure; + + if ( drawdown.value < 0 && wells().type[w] == INJECTOR) { + return false; + } + if ( drawdown.value > 0 && wells().type[w] == PRODUCER) { + return false; + } + } + return true; + } + void localInvert(Mat& istlA) const { for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) { @@ -484,7 +502,7 @@ namespace Opm { template void - computeWellFlux(const int& w, const double& Tw, const intensiveQuants& intQuants, const double& cdp, std::vector& cq_s) const + computeWellFlux(const int& w, const double& Tw, const intensiveQuants& intQuants, const double& cdp, const bool& allow_cf, std::vector& cq_s) const { const Opm::PhaseUsage& pu = fluid_->phaseUsage(); EvalWell bhp = getBhp(w); @@ -514,7 +532,7 @@ namespace Opm { if ( drawdown.value > 0 ) { //Do nothing if crossflow is not allowed - if (!wells().allow_cf[w] && wells().type[w] == INJECTOR) + if (!allow_cf && wells().type[w] == INJECTOR) return; // compute phase volumetric rates at standard conditions std::vector cq_ps(np, 0.0); @@ -539,7 +557,7 @@ namespace Opm { } else { //Do nothing if crossflow is not allowed - if (!wells().allow_cf[w] && wells().type[w] == PRODUCER) + if (!allow_cf && wells().type[w] == PRODUCER) return; // Using total mobilities