various fixes.

the most important one is that the formulation is changed in the recent
version
This commit is contained in:
Kai Bao
2019-11-27 13:05:11 +01:00
parent cc77c0e826
commit 3d7f0efe07
2 changed files with 32 additions and 18 deletions

View File

@@ -200,12 +200,13 @@ namespace mswellhelpers
// water in oil emulsion viscosity // water in oil emulsion viscosity
// TODO: maybe it should be two different ValueTypes. When we calculate the viscosity for transitional zone
template <typename ValueType> template <typename ValueType>
ValueType WIOEmulsionViscosity(const ValueType& oil_viscosity, const double water_liquid_fraction, ValueType WIOEmulsionViscosity(const ValueType& oil_viscosity, const ValueType& water_liquid_fraction,
const double max_visco_ratio) const double max_visco_ratio)
{ {
const double temp_value = 1. / (1. - (0.8415 / 0.7480 * water_liquid_fraction) ); const ValueType temp_value = 1. / (1. - (0.8415 / 0.7480 * water_liquid_fraction) );
const double viscosity_ratio = std::pow(temp_value, 2.5); const ValueType viscosity_ratio = Opm::pow(temp_value, 2.5);
if (viscosity_ratio <= max_visco_ratio) { if (viscosity_ratio <= max_visco_ratio) {
return oil_viscosity * viscosity_ratio; return oil_viscosity * viscosity_ratio;
@@ -220,11 +221,11 @@ namespace mswellhelpers
// oil in water emulsion viscosity // oil in water emulsion viscosity
template <typename ValueType> template <typename ValueType>
ValueType OIWEmulsionViscosity(const ValueType& water_viscosity, const double water_liquid_fraction, ValueType OIWEmulsionViscosity(const ValueType& water_viscosity, const ValueType& water_liquid_fraction,
const double max_visco_ratio) const double max_visco_ratio)
{ {
const double temp_value = 1. / (1. - (0.6019 / 0.6410) * (1. - water_liquid_fraction) ); const ValueType temp_value = 1. / (1. - (0.6019 / 0.6410) * (1. - water_liquid_fraction) );
const double viscosity_ratio = std::pow(temp_value, 2.5); const ValueType viscosity_ratio = Opm::pow(temp_value, 2.5);
if (viscosity_ratio <= max_visco_ratio) { if (viscosity_ratio <= max_visco_ratio) {
return water_viscosity * viscosity_ratio; return water_viscosity * viscosity_ratio;
@@ -239,8 +240,8 @@ namespace mswellhelpers
// calculating the viscosity of oil-water emulsion at local conditons // calculating the viscosity of oil-water emulsion at local conditons
template <typename ValueType> template <typename ValueType>
ValueType emulsionViscosity(const double water_fraction, const ValueType& water_viscosity, ValueType emulsionViscosity(const ValueType& water_fraction, const ValueType& water_viscosity,
const double oil_fraction, const ValueType& oil_viscosity, const ValueType& oil_fraction, const ValueType& oil_viscosity,
const SpiralICD& sicd) const SpiralICD& sicd)
{ {
const double width_transition = sicd.widthTransitionRegion(); const double width_transition = sicd.widthTransitionRegion();
@@ -251,16 +252,16 @@ namespace mswellhelpers
} }
const double critical_value = sicd.criticalValue(); const double critical_value = sicd.criticalValue();
const double transition_start_value = critical_value - width_transition / 2.0; const ValueType transition_start_value = critical_value - width_transition / 2.0;
const double transition_end_value = critical_value + width_transition / 2.0; const ValueType transition_end_value = critical_value + width_transition / 2.0;
const double liquid_fraction = water_fraction + oil_fraction; const ValueType liquid_fraction = water_fraction + oil_fraction;
// if there is no liquid, we just return zero // if there is no liquid, we just return zero
if (liquid_fraction == 0.) { if (liquid_fraction == 0.) {
return 0.; return 0.;
} }
const double water_liquid_fraction = water_fraction / liquid_fraction; const ValueType water_liquid_fraction = water_fraction / liquid_fraction;
const double max_visco_ratio = sicd.maxViscosityRatio(); const double max_visco_ratio = sicd.maxViscosityRatio();
if (water_liquid_fraction <= transition_start_value) { if (water_liquid_fraction <= transition_start_value) {

View File

@@ -3635,7 +3635,16 @@ namespace Opm
EvalWell pressure_equation = getSegmentPressure(seg); EvalWell pressure_equation = getSegmentPressure(seg);
pressure_equation = pressure_equation - pressureDropSpiralICD(seg); const int seg_upwind = upwinding_segments_[seg];
if (seg != seg_upwind) {
std::cout << " seg " << seg << " seg_upwind " << seg_upwind << std::endl;
}
const EvalWell pressure_drop_sicd = pressureDropSpiralICD(seg);
// std::cout << " pressure_drop_sicd for seg " << seg << " is " << pressure_drop_sicd.value()/1.e5 << std::endl;
// pressure_equation = pressure_equation - pressureDropSpiralICD(seg);
pressure_equation = pressure_equation - pressure_drop_sicd;
resWell_[seg][SPres] = pressure_equation.value(); resWell_[seg][SPres] = pressure_equation.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
@@ -3926,10 +3935,9 @@ namespace Opm
MultisegmentWell<TypeTag>:: MultisegmentWell<TypeTag>::
pressureDropSpiralICD(const int seg) const pressureDropSpiralICD(const int seg) const
{ {
// TODO: We have to consider the upwinding here
const SpiralICD& sicd = *segmentSet()[seg].spiralICD(); 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_fractions = segment_phase_fractions_[seg];
const std::vector<EvalWell>& phase_viscosities = segment_phase_viscosities_[seg]; const std::vector<EvalWell>& phase_viscosities = segment_phase_viscosities_[seg];
@@ -3961,8 +3969,8 @@ namespace Opm
// water_fraction + oil_fraction + gas_fraction should equal to one // water_fraction + oil_fraction + gas_fraction should equal to one
// calculating the water oil emulsion viscosity // calculating the water oil emulsion viscosity
// TODO: maybe we should keep the derivative of the fractions // TODO: maybe we should keep the derivative of the fractions
const EvalWell liquid_emulsion_viscosity = mswellhelpers::emulsionViscosity(water_fraction.value(), water_viscosity, const EvalWell liquid_emulsion_viscosity = mswellhelpers::emulsionViscosity(water_fraction, water_viscosity,
oil_fraction.value(), oil_viscosity, sicd); oil_fraction, oil_viscosity, sicd);
const EvalWell mixture_viscosity = (water_fraction + oil_fraction) * liquid_emulsion_viscosity + gas_fraction * gas_viscosities; 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 = segment_reservoir_volume_rates_[seg];
@@ -3974,10 +3982,15 @@ namespace Opm
using MathTool = MathToolbox<EvalWell>; using MathTool = MathToolbox<EvalWell>;
const EvalWell& density = segment_densities_[seg]; 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_value1 = MathTool::pow(density / density_cali, 0.75);
const EvalWell temp_value2 = MathTool::pow(mixture_viscosity / viscosity_cali, 0.25); 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; // const double base_strength = sicd.strength();// / density_cali;
// It looks like in 2016, they changed the formulation
const double strength = sicd.strength();
return temp_value1 * temp_value2 * strength * reservoir_rate_icd * reservoir_rate_icd;
} }