Use average density for hydro p-drop

This commit is contained in:
Stein Krogstad 2023-05-10 10:31:49 +02:00
parent 8ee3493aed
commit 0417f101bf
5 changed files with 23 additions and 13 deletions

View File

@ -219,14 +219,18 @@ assemblePressureLoss(const int seg,
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
assembleHydroPressureLoss(const int seg,
const int seg_upwind,
const EvalWell& hydro_pressure_drop,
const int seg_outlet,
const EvalWell& hydro_pressure_drop_seg,
const EvalWell& hydro_pressure_drop_outlet,
Equations& eqns1) const
{
MultisegmentWellEquationAccess<Scalar,numWellEq,Indices::numEq> eqns(eqns1);
eqns.residual()[seg][SPres] -= hydro_pressure_drop.value();
eqns.residual()[seg][SPres] -= (hydro_pressure_drop_seg.value() + hydro_pressure_drop_outlet.value());
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
eqns.D()[seg][seg][SPres][pv_idx] -= hydro_pressure_drop.derivative(pv_idx + Indices::numEq);
eqns.D()[seg][seg][SPres][pv_idx] -= hydro_pressure_drop_seg.derivative(pv_idx + Indices::numEq);
}
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
eqns.D()[seg][seg_outlet][SPres][pv_idx] -= hydro_pressure_drop_outlet.derivative(pv_idx + Indices::numEq);
}
}

View File

@ -89,9 +89,10 @@ public:
void assembleHydroPressureLoss(const int seg,
const int seg_upwind,
const EvalWell& hydro_pressure_drop,
Equations& eqs1) const;
const int seg_outlet,
const EvalWell& hydro_pressure_drop_seg,
const EvalWell& hydro_pressure_drop_outlet,
Equations& eqns1) const;
void assemblePressureEqExtraDerivatives(const int seg,
const int seg_upwind,

View File

@ -261,10 +261,13 @@ assembleDefaultPressureEq(const int seg,
}
// Do separately since different organization of derivatives
const auto hydro_pressure_drop = segments_.getHydroPressureLoss(seg);
const int seg_outlet = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const auto hydro_pressure_drop1 = segments_.getHalfHydroPressureLoss(seg, seg);
const auto hydro_pressure_drop2 = segments_.getHalfHydroPressureLoss(seg, seg_outlet);
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assembleHydroPressureLoss(seg, seg_upwind, hydro_pressure_drop, linSys_);
segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop.value();
assembleHydroPressureLoss(seg, seg_outlet, hydro_pressure_drop1, hydro_pressure_drop2, linSys_);
segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop1.value() + hydro_pressure_drop2.value();
}
template<typename FluidSystem, typename Indices, typename Scalar>

View File

@ -321,10 +321,11 @@ updateUpwindingSegments(const PrimaryVariables& primary_variables)
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
getHydroPressureLoss(const int seg) const
getHalfHydroPressureLoss(const int seg,
const int seg_side) const
{
//const int seg_upwind = upwinding_segments_[seg];
return densities_[seg] * well_.gravity() * depth_diffs_[seg];
return 0.5*densities_[seg_side] * well_.gravity() * depth_diffs_[seg];
}

View File

@ -51,7 +51,8 @@ public:
//! \brief Update upwinding segments.
void updateUpwindingSegments(const PrimaryVariables& primary_variables);
EvalWell getHydroPressureLoss(const int seg) const;
EvalWell getHalfHydroPressureLoss(const int seg,
const int seg_side) const;
//! Pressure difference between segment and perforation.
Scalar getPressureDiffSegPerf(const int seg,