OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

move pressureDropSpiralICD to MultisegmentWellSegments

Browse Source
This commit is contained in:
Arne Morten Kvarving 2022-12-19 14:20:55 +01:00
parent 4acfec9d1d
commit 133f2a92bb
4 changed files with 83 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
opm/simulators/wells/MultisegmentWellEval.cpp
View File

@ -179,86 +179,6 @@ extendEval(const Eval& in) const
return out;
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
pressureDropSpiralICD(const int seg) const
{
const SICD& sicd = this->segmentSet()[seg].spiralICD();
const int seg_upwind = segments_.upwinding_segments_[seg];
const std::vector<EvalWell>& phase_fractions = segments_.phase_fractions_[seg_upwind];
const std::vector<EvalWell>& phase_viscosities = segments_.phase_viscosities_[seg_upwind];
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_viscosity = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
gas_fraction = phase_fractions[gas_pos];
gas_viscosity = phase_viscosities[gas_pos];
}
EvalWell density = segments_.densities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
water_fraction.clearDerivatives();
water_viscosity.clearDerivatives();
oil_fraction.clearDerivatives();
oil_viscosity.clearDerivatives();
gas_fraction.clearDerivatives();
gas_viscosity.clearDerivatives();
density.clearDerivatives();
}
const EvalWell liquid_fraction = water_fraction + oil_fraction;
// viscosity contribution from the liquid
const EvalWell liquid_viscosity_fraction = liquid_fraction < 1.e-30 ? oil_fraction * oil_viscosity + water_fraction * water_viscosity :
liquid_fraction * mswellhelpers::emulsionViscosity(water_fraction, water_viscosity, oil_fraction, oil_viscosity, sicd);
const EvalWell mixture_viscosity = liquid_viscosity_fraction + gas_fraction * gas_viscosity;
const EvalWell reservoir_rate = segments_.mass_rates_[seg] / density;
const EvalWell reservoir_rate_icd = reservoir_rate * sicd.scalingFactor();
const double viscosity_cali = sicd.viscosityCalibration();
using MathTool = MathToolbox<EvalWell>;
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;
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
@ -499,7 +419,7 @@ assembleICDPressureEq(const int seg,
EvalWell icd_pressure_drop;
switch(this->segmentSet()[seg].segmentType()) {
case Segment::SegmentType::SICD :
icd_pressure_drop = pressureDropSpiralICD(seg);
icd_pressure_drop = segments_.pressureDropSpiralICD(seg);
break;
case Segment::SegmentType::AICD :
icd_pressure_drop = pressureDropAutoICD(seg, unit_system);

9
opm/simulators/wells/MultisegmentWellEval.hpp
View File

@ -36,10 +36,6 @@
#include <utility>
#include <vector>
namespace Dune {
template<class Matrix> class UMFPack;
}
namespace Opm
{
@ -129,9 +125,6 @@ protected:
EvalWell pressureDropAutoICD(const int seg,
const UnitSystem& unit_system) const;
// pressure drop for Spiral ICD segment (WSEGSICD)
EvalWell pressureDropSpiralICD(const int seg) const;
// pressure drop for sub-critical valve (WSEGVALV)
EvalWell pressureDropValve(const int seg) const;
@ -141,9 +134,7 @@ protected:
const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif_;
Equations linSys_; //!< The equation system
PrimaryVariables primary_variables_; //!< The primary variables
MSWSegments segments_; //!< Segment properties
// depth difference between perforations and the perforated grid cells

79
opm/simulators/wells/MultisegmentWellSegments.cpp
View File

@ -469,6 +469,85 @@ getFrictionPressureLoss(const int seg) const
return sign * mswellhelpers::frictionPressureLoss(length, diameter, area, roughness, density, mass_rate, visc);
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
pressureDropSpiralICD(const int seg) const
{
const auto& segment_set = well_.wellEcl().getSegments();
const SICD& sicd = segment_set[seg].spiralICD();
const int seg_upwind = upwinding_segments_[seg];
const std::vector<EvalWell>& phase_fractions = phase_fractions_[seg_upwind];
const std::vector<EvalWell>& phase_viscosities = phase_viscosities_[seg_upwind];
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_viscosity = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
gas_fraction = phase_fractions[gas_pos];
gas_viscosity = phase_viscosities[gas_pos];
}
EvalWell density = densities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
water_fraction.clearDerivatives();
water_viscosity.clearDerivatives();
oil_fraction.clearDerivatives();
oil_viscosity.clearDerivatives();
gas_fraction.clearDerivatives();
gas_viscosity.clearDerivatives();
density.clearDerivatives();
}
const EvalWell liquid_fraction = water_fraction + oil_fraction;
// viscosity contribution from the liquid
const EvalWell liquid_viscosity_fraction = liquid_fraction < 1.e-30 ? oil_fraction * oil_viscosity + water_fraction * water_viscosity :
liquid_fraction * mswellhelpers::emulsionViscosity(water_fraction, water_viscosity, oil_fraction, oil_viscosity, sicd);
const EvalWell mixture_viscosity = liquid_viscosity_fraction + gas_fraction * gas_viscosity;
const EvalWell reservoir_rate = mass_rates_[seg] / density;
const EvalWell reservoir_rate_icd = reservoir_rate * sicd.scalingFactor();
const double viscosity_cali = sicd.viscosityCalibration();
using MathTool = MathToolbox<EvalWell>;
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;
}
#define INSTANCE(...) \
template class MultisegmentWellSegments<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;

3
opm/simulators/wells/MultisegmentWellSegments.hpp
View File

@ -57,6 +57,9 @@ public:
EvalWell getFrictionPressureLoss(const int seg) const;
// pressure drop for Spiral ICD segment (WSEGSICD)
EvalWell pressureDropSpiralICD(const int seg) const;
// TODO: trying to use the information from the Well opm-parser as much
// as possible, it will possibly be re-implemented later for efficiency reason.