importing the old WSEGSICD implementation

with small adjustments to make it compile

opm/simulators/wells/MSWellHelpers.hpp

@@ -25,6 +25,7 @@
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/simulators/utils/DeferredLogger.hpp>
 #include <opm/common/ErrorMacros.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/MSW/SpiralICD.hpp>
 #include <dune/istl/solvers.hh>
 #if HAVE_UMFPACK
 #include <dune/istl/umfpack.hh>
@@ -197,6 +198,84 @@ namespace mswellhelpers
     }
 
 
+
+    // water in oil emulsion viscosity
+    template <typename ValueType>
+    ValueType WIOEmulsionViscosity(const ValueType& oil_viscosity, const double water_liquid_fraction,
+                                   const double max_visco_ratio)
+    {
+        const double temp_value = 1. / (1. - (0.8415 / 0.7480 * water_liquid_fraction) );
+        const double viscosity_ratio = std::pow(temp_value, 2.5);
+
+        if (viscosity_ratio <= max_visco_ratio) {
+            return oil_viscosity * viscosity_ratio;
+        } else {
+            return oil_viscosity * max_visco_ratio;
+        }
+    }
+
+
+
+
+
+    // oil in water emulsion viscosity
+    template <typename ValueType>
+    ValueType OIWEmulsionViscosity(const ValueType& water_viscosity, const double water_liquid_fraction,
+                                   const double max_visco_ratio)
+    {
+        const double temp_value = 1. / (1. - (0.6019 / 0.6410) * (1. - water_liquid_fraction) );
+        const double viscosity_ratio = std::pow(temp_value, 2.5);
+
+        if (viscosity_ratio <= max_visco_ratio) {
+            return water_viscosity * viscosity_ratio;
+        } else {
+            return water_viscosity * max_visco_ratio;
+        }
+    }
+
+
+
+
+
+    // calculating the viscosity of oil-water emulsion at local conditons
+    template <typename ValueType>
+    ValueType emulsionViscosity(const double water_fraction, const ValueType& water_viscosity,
+                                const double oil_fraction, const ValueType& oil_viscosity,
+                                const SpiralICD& sicd)
+    {
+        const double width_transition = sicd.widthTransitionRegion();
+
+        // it is just for now, we should be able to treat it.
+        if (width_transition <= 0.) {
+            OPM_THROW(std::runtime_error, "Not handling non-positive transition width now");
+        }
+
+        const double critical_value = sicd.criticalValue();
+        const double transition_start_value = critical_value - width_transition / 2.0;
+        const double transition_end_value = critical_value + width_transition / 2.0;
+
+        const double liquid_fraction = water_fraction + oil_fraction;
+        // if there is no liquid, we just return zero
+        if (liquid_fraction == 0.) {
+            return 0.;
+        }
+
+        const double water_liquid_fraction = water_fraction / liquid_fraction;
+
+        const double max_visco_ratio = sicd.maxViscosityRatio();
+        if (water_liquid_fraction <= transition_start_value) {
+            return WIOEmulsionViscosity(oil_viscosity, water_liquid_fraction, max_visco_ratio);
+        } else if(water_liquid_fraction >= transition_end_value) {
+            return OIWEmulsionViscosity(water_viscosity, water_liquid_fraction, max_visco_ratio);
+        } else { // in the transition region
+            const ValueType viscosity_start_transition = WIOEmulsionViscosity(oil_viscosity, transition_start_value, max_visco_ratio);
+            const ValueType viscosity_end_transition = OIWEmulsionViscosity(water_viscosity, transition_end_value, max_visco_ratio);
+            const ValueType emulsion_viscosity = (viscosity_start_transition * (transition_end_value - water_liquid_fraction)
+                                               + viscosity_end_transition * (water_liquid_fraction - transition_start_value) ) / width_transition;
+            return emulsion_viscosity;
+        }
+    }
+
 } // namespace mswellhelpers
 
 }

opm/simulators/wells/MultisegmentWell.hpp

@@ -273,6 +273,14 @@ namespace Opm
 
         mutable int debug_cost_counter_ = 0;
 
+        // TODO: this is the old implementation, it is possible the new value does not need it anymore
+        std::vector<EvalWell> segment_reservoir_volume_rates_;
+
+        std::vector<std::vector<EvalWell>> segment_phase_fractions_;
+
+        std::vector<std::vector<EvalWell>> segment_phase_viscosities_;
+
+
         void initMatrixAndVectors(const int num_cells) const;
 
         // protected functions
@@ -453,6 +461,7 @@ namespace Opm
         // be able to produce/inject .
         bool allDrawDownWrongDirection(const Simulator& ebos_simulator) const;
 
+
         boost::optional<double> computeBhpAtThpLimitProd(const Simulator& ebos_simulator,
                                                          const std::vector<Scalar>& B_avg,
                                                          const SummaryState& summary_state,
@@ -464,6 +473,13 @@ namespace Opm
                                                         DeferredLogger& deferred_logger) const;
 
         double maxPerfPress(const Simulator& ebos_simulator) const;
+
+        void assembleSICDPressureEq(const int seg) const;
+
+        void calculateFlowScalingFactors();
+
+        EvalWell pressureDropSpiralICD(const int seg) const;
+
     };
 
 }

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -49,6 +49,9 @@ namespace Opm
     , segment_mass_rates_(numberOfSegments(), 0.0)
     , segment_depth_diffs_(numberOfSegments(), 0.0)
     , upwinding_segments_(numberOfSegments(), 0)
+    , segment_reservoir_volume_rates_(numberOfSegments(), 0.0)
+    , segment_phase_fractions_(numberOfSegments(), std::vector<EvalWell>(num_components_, 0.0)) // number of phase here?
+    , segment_phase_viscosities_(numberOfSegments(), std::vector<EvalWell>(num_components_, 0.0)) // number of phase here?
     {
         // not handling solvent or polymer for now with multisegment well
         if (has_solvent) {
@@ -107,6 +110,10 @@ namespace Opm
             const double outlet_depth = outlet_segment.depth();
             segment_depth_diffs_[seg] = segment_depth - outlet_depth;
         }
+
+
+        // update the flow scaling factors for sicd segments
+        calculateFlowScalingFactors();
     }
 
 
@@ -1556,6 +1563,8 @@ namespace Opm
                 }
             }
 
+            segment_phase_viscosities_[seg] = visc;
+
             std::vector<EvalWell> mix(mix_s);
             if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                 const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
@@ -1581,12 +1590,22 @@ namespace Opm
                 segment_viscosities_[seg] += visc[comp_idx] * comp_fraction;
             }
 
+
+
             EvalWell density(0.0);
+            EvalWell surface_volume_rate(0.0);
+
             for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                 density += surf_dens[comp_idx] * mix_s[comp_idx];
+                surface_volume_rate += getSegmentRate(seg, comp_idx);
+
             }
             segment_densities_[seg] = density / volrat;
 
+            segment_reservoir_volume_rates_[seg] = surface_volume_rate / volrat;
+
+            segment_phase_fractions_[seg] = mix;
+
             // calculate the mass rates
             // TODO: for now, we are not considering the upwinding for this amount
             // since how to address the fact that the derivatives is not trivial for now
@@ -2884,10 +2903,19 @@ namespace Opm
                 const Opm::Schedule& schedule = ebosSimulator.vanguard().schedule();
                 assembleControlEq(well_state, schedule, summaryState, inj_controls, prod_controls, deferred_logger);
             } else {
-                assemblePressureEq(seg);
+                // TODO: maybe the following should go to the function assemblePressureEq()
+                if (segmentSet()[seg].segmentType() == Segment::SegmentType::SICD) {
+                    assembleSICDPressureEq(seg);
+                } else {
+                    // regular segment
+                    assemblePressureEq(seg);
+                }
             }
         }
     }
+
+
+
     template<typename TypeTag>
     bool
     MultisegmentWell<TypeTag>::
@@ -3586,11 +3614,48 @@ namespace Opm
             deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
                                     "Robust bhp(thp) solve failed for well " + name());
             return boost::optional<double>();
+	}
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    assembleSICDPressureEq(const int seg) const
+    {
+        // top segment can not be a spiral ICD device
+        assert(seg != 0);
+
+        // the pressure equation is something like
+        // p_seg - deltaP - p_outlet = 0.
+        // the major part is how to calculate the deltaP
+
+        EvalWell pressure_equation = getSegmentPressure(seg);
+
+        pressure_equation = pressure_equation - pressureDropSpiralICD(seg);
+
+        resWell_[seg][SPres] = pressure_equation.value();
+        for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+            duneD_[seg][seg][SPres][pv_idx] = pressure_equation.derivative(pv_idx + numEq);
+        }
+
+        // contribution from the outlet segment
+        const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
+        const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);
+
+        resWell_[seg][SPres] -= outlet_pressure.value();
+        for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+            duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq);
+
         }
     }
 
 
 
+
     template<typename TypeTag>
     boost::optional<double>
     MultisegmentWell<TypeTag>::
@@ -3792,6 +3857,48 @@ namespace Opm
 
 
 
+
+
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    calculateFlowScalingFactors()
+    {
+        // top segment will not be spiral ICD segment
+        for (int seg = 1; seg < numberOfSegments(); ++seg) {
+            const Segment& segment = segmentSet()[seg];
+            if (segment.segmentType() == Segment::SegmentType::SICD) {
+                SpiralICD& sicd = *segment.spiralICD();
+
+                // getting the segment length related to this ICD
+                const int parental_segment_number = segmentSet()[seg].outletSegment();
+                const double segment_length = segmentSet().segmentLength(parental_segment_number);
+
+                // getting the total completion length related to this ICD
+                // it should be connections
+                const auto& connections = well_ecl_.getConnections();
+                double total_connection_length = 0.;
+                for (const int conn : segment_perforations_[seg]) {
+                    const auto& connection = connections.get(conn);
+                    const double connection_length = connection.getSegDistEnd() - connection.getSegDistStart();
+                    assert(connection_length > 0.);
+                    total_connection_length += connection_length;
+                }
+
+                if (total_connection_length == 0.) {
+                    OPM_THROW(std::runtime_error, "zero total connection length is obtained when calcualting scaling factor");
+                }
+
+                sicd.updateScalingFactor(segment_length, total_connection_length);
+            }
+        }
+    }
+
+
+
+
+
+
     template<typename TypeTag>
     double
     MultisegmentWell<TypeTag>::
@@ -3812,4 +3919,66 @@ namespace Opm
     }
 
 
+
+
+
+    typename MultisegmentWell<TypeTag>::EvalWell
+    MultisegmentWell<TypeTag>::
+    pressureDropSpiralICD(const int seg) const
+    {
+        const SpiralICD& sicd = *segmentSet()[seg].spiralICD();
+
+        const double density_cali = sicd.densityCalibration();
+        const double base_strength = sicd.strength() / density_cali;
+
+        const std::vector<EvalWell>& phase_fractions = segment_phase_fractions_[seg];
+        const std::vector<EvalWell>& phase_viscosities = segment_phase_viscosities_[seg];
+
+        EvalWell water_fraction = 0.;
+        EvalWell water_viscosity = 0.;
+        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+            const int water_pos = phaseUsage().phase_pos[Water];
+            water_fraction = phase_fractions[water_pos];
+            water_viscosity = phase_viscosities[water_pos];
+        }
+
+        EvalWell oil_fraction = 0.;
+        EvalWell oil_viscosity = 0.;
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+            const int oil_pos = phaseUsage().phase_pos[Oil];
+            oil_fraction = phase_fractions[oil_pos];
+            oil_viscosity = phase_viscosities[oil_pos];
+        }
+
+        EvalWell gas_fraction = 0.;
+        EvalWell gas_viscosities = 0.;
+        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            const int gas_pos = phaseUsage().phase_pos[Gas];
+            gas_fraction = phase_fractions[gas_pos];
+            gas_viscosities = phase_viscosities[gas_pos];
+        }
+
+        // water_fraction + oil_fraction + gas_fraction should equal to one
+        // calculating the water oil emulsion viscosity
+        // TODO: maybe we should keep the derivative of the fractions
+        const EvalWell liquid_emulsion_viscosity = mswellhelpers::emulsionViscosity(water_fraction.value(), water_viscosity,
+                                                 oil_fraction.value(), oil_viscosity, sicd);
+        const EvalWell mixture_viscosity = (water_fraction + oil_fraction) * liquid_emulsion_viscosity + gas_fraction * gas_viscosities;
+
+        const EvalWell& reservoir_rate = segment_reservoir_volume_rates_[seg];
+
+        const EvalWell reservoir_rate_icd = reservoir_rate * sicd.scalingFactor();
+
+        const double viscosity_cali = sicd.viscosityCalibration();
+
+        using MathTool = MathToolbox<EvalWell>;
+
+        const EvalWell& density = segment_densities_[seg];
+        const EvalWell temp_value1 = MathTool::pow(density / density_cali, 0.75);
+        const EvalWell temp_value2 = MathTool::pow(mixture_viscosity / viscosity_cali, 0.25);
+
+        return temp_value1 * temp_value2 * base_strength * reservoir_rate_icd * reservoir_rate_icd;
+    }
+
+
 }