mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge pull request #2178 from GitPaean/support_wsegsicd_rebased
Support WSEGSICD implementation
This commit is contained in:
Atgeirr Flø Rasmussen
GitHub
commit
4bad782a1f
@ -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
|
||||
|
||||
@ -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,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 <typename ValueType>
|
||||
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 <typename ValueType>
|
||||
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 <typename ValueType>
|
||||
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
|
||||
|
||||
}
|
||||
|
||||
@ -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,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;
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
@ -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,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<EvalWell> 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 {
|
||||
// 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>::
|
||||
@ -3591,6 +3614,44 @@ namespace Opm
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename TypeTag>
|
||||
void
|
||||
MultisegmentWell<TypeTag>::
|
||||
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<typename TypeTag>
|
||||
boost::optional<double>
|
||||
MultisegmentWell<TypeTag>::
|
||||
@ -3792,6 +3853,43 @@ namespace Opm
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename TypeTag>
|
||||
void
|
||||
MultisegmentWell<TypeTag>::
|
||||
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<typename TypeTag>
|
||||
double
|
||||
MultisegmentWell<TypeTag>::
|
||||
@ -3812,4 +3910,70 @@ namespace Opm
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename TypeTag>
|
||||
typename MultisegmentWell<TypeTag>::EvalWell
|
||||
MultisegmentWell<TypeTag>::
|
||||
pressureDropSpiralICD(const int seg) const
|
||||
{
|
||||
// TODO: We have to consider the upwinding here
|
||||
const SpiralICD& sicd = *segmentSet()[seg].spiralICD();
|
||||
|
||||
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 = 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<EvalWell>;
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
@ -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=""
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user