handling the acceleration pressure drop

not tested yet.

opm/autodiff/MSWellHelpers.hpp

@@ -119,9 +119,9 @@ namespace mswellhelpers
     // density is density
     // roughness is the absolute roughness
     // mu is the average phase viscosity
-    template <class ValueType>
-    ValueType frictionPressureLoss(const double l, const double diameter, const double area, const ValueType& density,
-                                   const ValueType& w, const double roughness, const ValueType& mu)
+    template <typename ValueType>
+    ValueType frictionPressureLoss(const double l, const double diameter, const double area, const double roughness,
+                                   const ValueType& density, const ValueType& w, const ValueType& mu)
     {
         const double f = calculateFrictionFactor(area, diameter, w.value(), roughness, mu.value());
         return f * l * w * w / (area * area * diameter * density);
@@ -130,6 +130,14 @@ namespace mswellhelpers
 
 
 
+
+    template <typename ValueType>
+    ValueType velocityHead(const double area, const ValueType& mass_rate, const ValueType& density)
+    {
+        return (0.5 * mass_rate * mass_rate / (area * area * density));
+    }
+
+
 } // namespace mswellhelpers
 
 }

opm/autodiff/MultisegmentWell.hpp

@@ -335,12 +335,18 @@ namespace Opm
         // frictinal pressure loss
         EvalWell getFrictionPressureLoss(const int seg) const;
 
+        void handleAccelerationPressureLoss(const int seg) const;
+
         // handling the overshooting and undershooting of the fractions
         void processFractions(const int seg) const;
 
         void updateWellStateFromPrimaryVariables(WellState& well_state) const;
 
         double scalingFactor(const int comp_idx) const;
+
+        bool frictionalPressureLossConsidered() const;
+
+        bool accelerationalPressureLossConsidered() const;
     };
 
 }

opm/autodiff/MultisegmentWell_impl.hpp

@@ -1524,6 +1524,11 @@ namespace Opm
         // we only consider the hydrostatic pressure loss first
         pressure_equation -= getHydroPressureLoss(seg);
 
+        if (frictionalPressureLossConsidered()) {
+            // TODO: deciding the direction of friction later
+            pressure_equation -= getFrictionPressureLoss(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);
@@ -1537,6 +1542,10 @@ namespace Opm
         for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
             duneD_[seg][outlet_segment_location][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq);
         }
+
+        if (accelerationalPressureLossConsidered()) {
+            handleAccelerationPressureLoss(seg);
+        }
     }
 
 
@@ -1570,7 +1579,49 @@ namespace Opm
         const double area = segmentSet()[seg].crossArea();
         const double diameter = segmentSet()[seg].internalDiameter();
 
-        return frictionPressureLoss(length, diameter, area, density, mass_rate, roughness, visc);
+        return mswellhelpers::frictionPressureLoss(length, diameter, area, roughness, density, mass_rate, visc);
+    }
+
+
+
+
+
+    template <typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    handleAccelerationPressureLoss(const int seg) const
+    {
+        // handle the out velcocity head
+        const double area = segmentSet()[seg].crossArea();
+        const EvalWell mass_rate = segment_mass_rates_[seg];
+        const EvalWell density = segment_densities_[seg];
+        const EvalWell out_velocity_head = mswellhelpers::velocityHead(area, mass_rate, density);
+
+        // TODO: the sign is really hard and not sure
+        resWell_[seg][SPres] -= out_velocity_head.value();
+        for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+            duneD_[seg][seg][SPres][pv_idx] -= out_velocity_head.derivative(pv_idx + numEq);
+        }
+
+        // calcuate the maximum cross-area among the segment and its inlet segments
+        double max_area = area;
+        for (const int inlet : segment_inlets_[seg]) {
+            const double inlet_area = segmentSet()[seg].crossArea();
+            if (inlet_area > max_area) {
+                max_area = inlet_area;
+            }
+        }
+
+        // handling the velocity head of intlet segments
+        for (const int inlet : segment_inlets_[seg]) {
+            const EvalWell density = segment_densities_[inlet];
+            const EvalWell mass_rate = segment_mass_rates_[inlet];
+            const EvalWell inlet_velocity_head = mswellhelpers::velocityHead(area, mass_rate, density);
+            resWell_[seg][SPres] += inlet_velocity_head.value();
+            for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+                duneD_[seg][inlet][SPres][pv_idx] += inlet_velocity_head.derivative(pv_idx + numEq);
+            }
+        }
     }
 
 
@@ -1711,7 +1762,8 @@ namespace Opm
 
     template <typename TypeTag>
     double
-    MultisegmentWell<TypeTag>::scalingFactor(const int comp_idx) const
+    MultisegmentWell<TypeTag>::
+    scalingFactor(const int comp_idx) const
     {
         const double* distr = well_controls_get_current_distr(well_controls_);
 
@@ -1736,4 +1788,30 @@ namespace Opm
         assert(false);
         return 1.0;
     }
+
+
+
+
+
+    template <typename TypeTag>
+    bool
+    MultisegmentWell<TypeTag>::
+    frictionalPressureLossConsidered() const
+    {
+        // HF- and HFA needs to consider frictional pressure loss
+        return (SegmentSet().compPressureDrop() != WellSegment::H__);
+    }
+
+
+
+
+
+    template <typename TypeTag>
+    bool
+    MultisegmentWell<TypeTag>::
+    accelerationalPressureLossConsidered() const
+    {
+        return (SegmentSet().compPressureDrop() == WellSegment::HFA);
+    }
+
 }