From 61ed33c4581e71f3ecd81bf71db26c8dd7de66f0 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 29 Nov 2019 17:23:26 +0100 Subject: [PATCH] implementing the support for WSEGVALV --- opm/simulators/wells/MSWellHelpers.hpp | 10 +- opm/simulators/wells/MultisegmentWell.hpp | 4 + .../wells/MultisegmentWell_impl.hpp | 172 +++++++++++++----- 3 files changed, 138 insertions(+), 48 deletions(-) diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp index d070dd7b0..811906f31 100644 --- a/opm/simulators/wells/MSWellHelpers.hpp +++ b/opm/simulators/wells/MSWellHelpers.hpp @@ -188,7 +188,15 @@ namespace mswellhelpers } - + template + ValueType valveContrictionPressureLoss(const ValueType& mass_rate, const ValueType& density, + const double area_con, const double cv) + { + // the formulation is adjusted a little bit for convinience + // velocity = mass_rate / (density * area) is applied to the original formulation + const double area = (area_con > 1.e-10 ? area_con : 1.e-10); + return mass_rate * mass_rate / (2. * density * cv * cv * area * area); + } template diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index cb256475f..8ae3a94a8 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -481,6 +481,10 @@ namespace Opm EvalWell pressureDropSpiralICD(const int seg) const; + // assemble the pressure equation for sub-critical valve (WSEGVALV) + void assembleValvePressureEq(const int seg) const; + + EvalWell pressureDropValve(const int seg) const; }; } diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index aebfef64a..1c0860a04 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -21,6 +21,7 @@ #include #include +#include namespace Opm { @@ -2654,7 +2655,7 @@ namespace Opm if (seg == 0) { // top segment well_state.bhp()[index_of_well_] = well_state.segPress()[seg + top_segment_index]; } - } + } updateThp(well_state, deferred_logger); } @@ -2754,7 +2755,7 @@ namespace Opm std::ostringstream sstr; if (is_stagnate) { sstr << " well " << name() << " observes stagnation in inner iteration " << it << "\n"; - + } if (is_oscillate) { sstr << " well " << name() << " observes oscillation in inner iteration " << it << "\n"; @@ -2936,12 +2937,16 @@ namespace Opm 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); - } + switch(segmentSet()[seg].segmentType()) { + case Segment::SegmentType::SICD : + assembleSICDPressureEq(seg); + break; + case Segment::SegmentType::VALVE : + assembleValvePressureEq(seg); + break; + default : + assemblePressureEq(seg); + } } } } @@ -3452,6 +3457,88 @@ namespace Opm + + + + + + + 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 + void + MultisegmentWell:: + assembleValvePressureEq(const int seg) const + { + // TODO: upwinding needs to be taken care of + // top segment can not be a spiral ICD device + assert(seg != 0); + + // const Valve& valve = *segmentSet()[seg].Valve(); + + // 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); + + // const int seg_upwind = upwinding_segments_[seg]; + + pressure_equation = pressure_equation - pressureDropValve(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:: @@ -3652,43 +3739,6 @@ namespace Opm - - 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:: @@ -3925,8 +3975,6 @@ namespace Opm - - template double MultisegmentWell:: @@ -4013,4 +4061,34 @@ namespace Opm } + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + pressureDropValve(const int seg) const + { + const Valve& valve = *segmentSet()[seg].valve(); + + const EvalWell& mass_rate = segment_mass_rates_[seg]; + const EvalWell& visc = segment_viscosities_[seg]; + const EvalWell& density = segment_densities_[seg]; + const double additional_length = valve.pipeAdditionalLength(); + const double roughness = valve.pipeRoughness(); + const double diameter = valve.pipeDiameter(); + const double area = valve.pipeCrossArea(); + + const EvalWell friction_pressure_loss = + mswellhelpers::frictionPressureLoss(additional_length, diameter, area, roughness, density, mass_rate, visc); + + const double area_con = valve.conCrossArea(); + const double cv = valve.conFlowCoefficient(); + + const EvalWell constriction_pressure_loss = + mswellhelpers::valveContrictionPressureLoss(mass_rate, density, area_con, cv); + + const double sign = mass_rate <= 0. ? 1.0 : -1.0; + return sign * (friction_pressure_loss + constriction_pressure_loss); + } + }