Merge pull request #4346 from GitPaean/support_cake_filtration

A simple cake model to simulate formation damage due to suspended solids in injection water

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -1421,6 +1421,21 @@ void BlackoilWellModelGeneric::initInjMult() {
     }
 }
 
+
+void BlackoilWellModelGeneric::updateFiltrationParticleVolume(const double dt, const size_t water_index) {
+    for (auto& well : this->well_container_generic_) {
+        if (well->isInjector() && well->wellEcl().getFilterConc() > 0.) {
+            auto &values =  this->filtration_particle_volume_[well->name()];
+            const auto& ws = this->wellState().well(well->indexOfWell());
+            if (values.empty()) {
+                values.assign(ws.perf_data.size(), 0.); // initializing to be zero
+            }
+            well->updateFiltrationParticleVolume(dt, water_index, this->wellState(), values);
+        }
+    }
+
+}
+
 void BlackoilWellModelGeneric::updateInjMult(DeferredLogger& deferred_logger) {
     for (const auto& well : this->well_container_generic_) {
         if (well->isInjector() && well->wellEcl().getInjMultMode() != Well::InjMultMode::NONE) {

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -376,6 +376,8 @@ protected:
 
     void updateInjMult(DeferredLogger& deferred_logger);
 
+    void updateFiltrationParticleVolume(const double dt, const size_t water_index);
+
     // create the well container
     virtual void createWellContainer(const int time_step) = 0;
     virtual void initWellContainer(const int reportStepIdx) = 0;
@@ -438,6 +440,9 @@ protected:
     // previous injection multiplier, it is used in the injection multiplier calculation for WINJMULT keyword
     std::unordered_map<std::string, std::vector<double>> prev_inj_multipliers_;
 
+    // the volume of the filtration particles during water injection
+    std::unordered_map<std::string, std::vector<double>> filtration_particle_volume_;
+
     /*
       The various wellState members should be accessed and modified
       through the accessor functions wellState(), prevWellState(),

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -310,6 +310,16 @@ namespace Opm {
             well->setGuideRate(&guideRate_);
         }
 
+        for (auto& well : well_container_) {
+            if (well->isInjector()) {
+                const auto it = this->filtration_particle_volume_.find(well->name());
+                if (it != this->filtration_particle_volume_.end()) {
+                    const auto& filtration_particle_volume = it->second;
+                    well->updateInjFCMult(filtration_particle_volume, local_deferredLogger);
+                }
+            }
+        }
+
         // Close completions due to economic reasons
         for (auto& well : well_container_) {
             well->closeCompletions(wellTestState());
@@ -474,6 +484,10 @@ namespace Opm {
             }
         }
 
+        if (Indices::waterEnabled) {
+            this->updateFiltrationParticleVolume(dt, FluidSystem::waterPhaseIdx);
+        }
+
         // at the end of the time step, updating the inj_multiplier saved in WellState for later use
         this->updateInjMult(local_deferredLogger);
 

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -24,6 +24,7 @@
 
 #include <opm/common/ErrorMacros.hpp>
 
+#include <opm/input/eclipse/Schedule/Well/FilterCake.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellBrineProperties.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellFoamProperties.hpp>
@@ -97,6 +98,9 @@ WellInterfaceGeneric::WellInterfaceGeneric(const Well& well,
             saturation_table_number_[perf] = pd.satnum_id;
             ++perf;
         }
+        if (this->isInjector()) {
+            inj_fc_multiplier_.resize(number_of_perforations_, 1.0);
+        }
     }
 
     // initialization of the completions mapping
@@ -714,4 +718,82 @@ checkNegativeWellPotentials(std::vector<double>& well_potentials,
     }
 }
 
+void WellInterfaceGeneric::
+updateFiltrationParticleVolume(const double dt, const size_t water_index,
+                               const WellState& well_state, std::vector<double>& filtration_particle_volume) const
+{
+    if (!this->isInjector()) {
+        return;
+    }
+
+    const auto injectorType = this->well_ecl_.injectorType();
+    if (injectorType != InjectorType::WATER) {
+        return;
+    }
+
+    const double conc = this->well_ecl_.getFilterConc();
+    if (conc == 0.) {
+        return;
+    }
+
+    // it is currently used for the filter cake modeling related to formation damage study
+    auto& ws = well_state.well(this->index_of_well_);
+    const auto& connection_rates = ws.perf_data.phase_rates;
+
+    const std::size_t np = well_state.numPhases();
+    for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
+        // not considering the production water
+        const double water_rates = std::max(0., connection_rates[perf * np + water_index]);
+        filtration_particle_volume[perf] += water_rates * conc * dt;
+    }
+}
+
+void WellInterfaceGeneric::
+updateInjFCMult(const std::vector<double>& filtration_particle_volume, DeferredLogger& deferred_logger) {
+    for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
+        const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
+        const auto& connections = this->well_ecl_.getConnections();
+        const auto& connection = connections[perf_ecl_index];
+        if (this->isInjector() && connection.filterCakeActive()) {
+            const auto& filter_cake = connection.getFilterCake();
+            const double area = connection.getFilterCakeArea();
+            const double poro = filter_cake.poro;
+            const double perm = filter_cake.perm;
+            const double rw = connection.getFilterCakeRadius();
+            const auto cr0 = connection.r0();
+            const auto crw = connection.rw();
+            const double K = connection.Kh() / connection.connectionLength();
+            const double factor = filter_cake.sf_multiplier;
+            // the thickness of the filtration cake
+            const double thickness = filtration_particle_volume[perf] / (area * (1. - poro));
+
+            double skin_factor = 0.;
+            switch (filter_cake.geometry) {
+                case FilterCake::FilterCakeGeometry::LINEAR: {
+                    skin_factor = thickness / rw * K / perm * factor;
+                    break;
+                }
+                case FilterCake::FilterCakeGeometry::RADIAL: {
+                    const double rc = std::sqrt(rw * rw + 2. * rw * thickness);
+                    skin_factor = K / perm * std::log(rc / rw) * factor;
+                    break;
+                }
+                default:
+                    OPM_DEFLOG_THROW(std::runtime_error,
+                                     fmt::format(" Invalid filtration cake geometry type ({}) for well {}",
+                                                 FilterCake::filterCakeGeometryToString(filter_cake.geometry), name()),
+                                     deferred_logger);
+            }
+            // the original skin factor for the connection
+            const auto cskinfactor = connection.skinFactor();
+            // compute a multiplier for the well connection transmissibility
+            const auto denom = std::log(cr0 / std::min(crw, cr0)) + cskinfactor;
+            const auto denom2 = denom + skin_factor;
+            this->inj_fc_multiplier_[perf] = denom / denom2;
+        } else {
+            this->inj_fc_multiplier_[perf] = 1.0;
+        }
+    }
+}
+
 } // namespace Opm

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -188,6 +188,13 @@ public:
     // it might change in the future
     double getInjMult(const int perf, const double bhp, const double perf_pres) const;
 
+    // update the water injection volume, it will be used for calculation related to cake filtration due to injection activity
+    void updateFiltrationParticleVolume(const double dt, const size_t water_index,
+                                        const WellState& well_state, std::vector<double>& filtration_particle_volume) const;
+
+    // update the multiplier for well transmissbility due to cake filteration
+    void updateInjFCMult(const std::vector<double>& filtration_particle_volume, DeferredLogger& deferred_logger);
+
     // whether a well is specified with a non-zero and valid VFP table number
     bool isVFPActive(DeferredLogger& deferred_logger) const;
 
@@ -369,6 +376,9 @@ protected:
     // which intends to keep the fracturing open
     std::vector<double> prev_inj_multiplier_;
 
+    // the multiplier due to injection filtration cake
+    std::vector<double> inj_fc_multiplier_;
+
     double well_efficiency_factor_;
     const VFPProperties* vfp_properties_;
     const GuideRate* guide_rate_;

opm/simulators/wells/WellInterface_impl.hpp

@@ -1267,6 +1267,17 @@ namespace Opm
                 OPM_DEFLOG_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent", deferred_logger);
             }
         }
+
+        if (this->isInjector()) {
+            const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
+            const auto& connections = this->well_ecl_.getConnections();
+            const auto& connection = connections[perf_ecl_index];
+            if (connection.filterCakeActive()) {
+                for (auto& val : mob) {
+                    val *= this->inj_fc_multiplier_[perf];
+                }
+            }
+        }
     }
 
 

parallelTests.cmake

@@ -181,6 +181,14 @@ add_test_compare_parallel_simulation(CASENAME winjmult_msw
                                      DIR winjmult
                                      TEST_ARGS --enable-tuning=true)
 
+add_test_compare_parallel_simulation(CASENAME winjdam_msw
+                                     FILENAME WINJDAM_MSW
+                                     SIMULATOR flow
+                                     ABS_TOL ${abs_tol}
+                                     REL_TOL 0.1
+                                     DIR winjdam
+                                     TEST_ARGS --enable-tuning=true)
+
 add_test_compare_parallel_simulation(CASENAME 3_a_mpi_multflt_mod2
                                      FILENAME 3_A_MPI_MULTFLT_SCHED_MODEL2
                                      SIMULATOR flow

regressionTests.cmake

@@ -1213,6 +1213,20 @@ add_test_compareECLFiles(CASENAME winjmult_msw
                          REL_TOL ${rel_tol}
                          DIR winjmult
                          TEST_ARGS --enable-tuning=true)
+add_test_compareECLFiles(CASENAME winjdam_stdw
+                         FILENAME WINJDAM_STDW
+                         SIMULATOR flow
+                         ABS_TOL ${abs_tol}
+                         REL_TOL ${rel_tol}
+                         DIR winjdam
+                         TEST_ARGS --enable-tuning=true)
+add_test_compareECLFiles(CASENAME winjdam_msw
+                         FILENAME WINJDAM_MSW
+                         SIMULATOR flow
+                         ABS_TOL ${abs_tol}
+                         REL_TOL ${rel_tol}
+                         DIR winjdam
+                         TEST_ARGS --enable-tuning=true)
 add_test_compareECLFiles(CASENAME 01_vappars
                          FILENAME VAPPARS-01
                          SIMULATOR flow