diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake index 9f8d04b42..5c28cff1b 100755 --- a/compareECLFiles.cmake +++ b/compareECLFiles.cmake @@ -382,6 +382,12 @@ add_test_compareECLFiles(CASENAME editnnc_model2 REL_TOL ${rel_tol} DIR model2) +add_test_compareECLFiles(CASENAME wsegsicd + FILENAME TEST_WSEGSICD + SIMULATOR flow + ABS_TOL ${abs_tol} + REL_TOL ${rel_tol}) + add_test_compareECLFiles(CASENAME nnc FILENAME NNC_AND_EDITNNC SIMULATOR flow diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp index 09851c289..d070dd7b0 100644 --- a/opm/simulators/wells/MSWellHelpers.hpp +++ b/opm/simulators/wells/MSWellHelpers.hpp @@ -25,6 +25,7 @@ #include #include #include +#include #include #if HAVE_UMFPACK #include @@ -197,6 +198,85 @@ namespace mswellhelpers } + + // water in oil emulsion viscosity + // TODO: maybe it should be two different ValueTypes. When we calculate the viscosity for transitional zone + template + ValueType WIOEmulsionViscosity(const ValueType& oil_viscosity, const ValueType& water_liquid_fraction, + const double max_visco_ratio) + { + const ValueType temp_value = 1. / (1. - (0.8415 / 0.7480 * water_liquid_fraction) ); + const ValueType viscosity_ratio = Opm::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 + ValueType OIWEmulsionViscosity(const ValueType& water_viscosity, const ValueType& water_liquid_fraction, + const double max_visco_ratio) + { + const ValueType temp_value = 1. / (1. - (0.6019 / 0.6410) * (1. - water_liquid_fraction) ); + const ValueType viscosity_ratio = Opm::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 + ValueType emulsionViscosity(const ValueType& water_fraction, const ValueType& water_viscosity, + const ValueType& 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 ValueType transition_start_value = critical_value - width_transition / 2.0; + const ValueType transition_end_value = critical_value + width_transition / 2.0; + + const ValueType liquid_fraction = water_fraction + oil_fraction; + // if there is no liquid, we just return zero + if (liquid_fraction == 0.) { + return 0.; + } + + const ValueType 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 } diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 6760c0b78..cb256475f 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/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 segment_reservoir_volume_rates_; + + std::vector> segment_phase_fractions_; + + std::vector> 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 computeBhpAtThpLimitProd(const Simulator& ebos_simulator, const std::vector& B_avg, const SummaryState& summary_state, @@ -464,6 +473,14 @@ namespace Opm DeferredLogger& deferred_logger) const; double maxPerfPress(const Simulator& ebos_simulator) const; + + void assembleSICDPressureEq(const int seg) const; + + // TODO: when more ICD devices join, we should have a better interface to do this + void calculateSICDFlowScalingFactors(); + + EvalWell pressureDropSpiralICD(const int seg) const; + }; } diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 193093415..50b24ac8b 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/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(num_components_, 0.0)) // number of phase here? + , segment_phase_viscosities_(numberOfSegments(), std::vector(num_components_, 0.0)) // number of phase here? { // not handling solvent or polymer for now with multisegment well if (has_solvent) { @@ -107,6 +110,9 @@ 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 + calculateSICDFlowScalingFactors(); } @@ -1556,16 +1562,20 @@ namespace Opm } } + segment_phase_viscosities_[seg] = visc; + std::vector mix(mix_s); if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const EvalWell d = 1.0 - rs * rv; + if (rs != 0.0) { // rs > 0.0? - mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / (1. - rs * rv); + mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d; } if (rv != 0.0) { // rv > 0.0? - mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / (1. - rs * rv); + mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d; } } @@ -1577,8 +1587,10 @@ namespace Opm segment_viscosities_[seg] = 0.; // calculate the average viscosity for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - const EvalWell comp_fraction = mix[comp_idx] / b[comp_idx] / volrat; - segment_viscosities_[seg] += visc[comp_idx] * comp_fraction; + const EvalWell fraction = mix[comp_idx] / b[comp_idx] / volrat; + // TODO: a little more work needs to be done to handle the negative fractions here + segment_phase_fractions_[seg][comp_idx] = fraction; // >= 0.0 ? fraction : 0.0; + segment_viscosities_[seg] += visc[comp_idx] * segment_phase_fractions_[seg][comp_idx]; } EvalWell density(0.0); @@ -1597,6 +1609,8 @@ namespace Opm const EvalWell rate = getSegmentRate(seg, comp_idx); segment_mass_rates_[seg] += rate * surf_dens[comp_idx]; } + + segment_reservoir_volume_rates_[seg] = segment_mass_rates_[seg] / segment_densities_[seg]; } } @@ -2884,10 +2898,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 bool MultisegmentWell:: @@ -3586,11 +3609,49 @@ namespace Opm deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE", "Robust bhp(thp) solve failed for well " + name()); return boost::optional(); + } + } + + + + + + template + void + MultisegmentWell:: + assembleSICDPressureEq(const int seg) const + { + // TODO: upwinding needs to be taken care of + // 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 boost::optional MultisegmentWell:: @@ -3792,6 +3853,43 @@ namespace Opm + + + template + void + MultisegmentWell:: + calculateSICDFlowScalingFactors() + { + // 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) { + // 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; + } + + SpiralICD& sicd = *segment.spiralICD(); + sicd.updateScalingFactor(segment_length, total_connection_length); + } + } + } + + + + + + template double MultisegmentWell:: @@ -3812,4 +3910,70 @@ namespace Opm } + + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + pressureDropSpiralICD(const int seg) const + { + // TODO: We have to consider the upwinding here + const SpiralICD& sicd = *segmentSet()[seg].spiralICD(); + + const std::vector& phase_fractions = segment_phase_fractions_[seg]; + const std::vector& phase_viscosities = segment_phase_viscosities_[seg]; + + EvalWell water_fraction = 0.; + EvalWell water_viscosity = 0.; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int water_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + 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 = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + 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 = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + gas_fraction = phase_fractions[gas_pos]; + gas_viscosities = phase_viscosities[gas_pos]; + } + + const EvalWell liquid_emulsion_viscosity = mswellhelpers::emulsionViscosity(water_fraction, water_viscosity, + oil_fraction, 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; + + const EvalWell& density = segment_densities_[seg]; + const double density_cali = sicd.densityCalibration(); + const EvalWell temp_value1 = MathTool::pow(density / density_cali, 0.75); + const EvalWell temp_value2 = MathTool::pow(mixture_viscosity / viscosity_cali, 0.25); + + // formulation before 2016, base_strength is used + // const double base_strength = sicd.strength() / density_cali; + // formulation since 2016, strength is used instead + const double strength = sicd.strength(); + + const double sign = reservoir_rate_icd <= 0. ? 1.0 : -1.0; + + return sign * temp_value1 * temp_value2 * strength * reservoir_rate_icd * reservoir_rate_icd; + } + + } diff --git a/tests/update_reference_data.sh b/tests/update_reference_data.sh index 6a96da64f..b8b92a30c 100755 --- a/tests/update_reference_data.sh +++ b/tests/update_reference_data.sh @@ -88,6 +88,7 @@ tests[polymer_injectivity]="flow polymer_injectivity 2D_POLYMER_INJECTIVITY" tests[nnc]="flow editnnc NNC_AND_EDITNNC nnc" tests[udq]="flow udq_actionx UDQ_WCONPROD udq_wconprod" tests[spe1_foam]="flow spe1_foam SPE1FOAM" +tests[wsegsicd]="flow wsegsicd TEST_WSEGSICD" changed_tests=""