Merge branch 'fix/navier-stokes-rot-sym' into 'master'

Fix/navier stokes rot sym

Closes #927

See merge request !2236

dumux/discretization/extrusion.hh

@@ -53,10 +53,13 @@ struct NoExtrusion
 /*!
  * \ingroup Discretization
  * \brief Rotation symmetric extrusion policy for rotating about an external axis
+ * \tparam radAx The radial axis perpendicular to the symmetry axis (0 = x, 1 = y)
  */
-template<int radialAxis = 0>
+template<int radAx = 0>
 struct RotationalExtrusion
 {
+    static constexpr int radialAxis = radAx;
+
     /*!
      * \brief Transformed sub-control-volume face area
      * \note Mid-point rule integrals are only exact for constants

dumux/freeflow/navierstokes/staggered/localresidual.hh

@@ -169,6 +169,24 @@ public:
         source += problem.gravity()[scvf.directionIndex()] * insideVolVars.density();
         source += problem.source(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())];
 
+        // Axisymmetric problems in 2D feature an extra source terms arising from the transformation to cylindrical coordinates.
+        // See Ferziger/Peric: Computational methods for fluid dynamics chapter 8.
+        // https://doi.org/10.1007/978-3-540-68228-8 (page 301)
+        if constexpr (ModelTraits::dim() == 2 && std::is_same_v<Extrusion, RotationalExtrusion<Extrusion::radialAxis>>)
+        {
+            if (scvf.directionIndex() == Extrusion::radialAxis)
+            {
+                // Velocity term
+                const auto r = scvf.center()[scvf.directionIndex()];
+                source -= -2.0*insideVolVars.effectiveViscosity() * elemFaceVars[scvf].velocitySelf() / (r*r);
+
+                // Pressure term (needed because we incorporate pressure in terms of a surface integral).
+                // grad(p) becomes div(pI) + (p/r)*n_r in cylindrical coordinates. The second term
+                // is new with respect to Cartesian coordinates and handled below as a source term.
+                source += insideVolVars.pressure()/r;
+            }
+        }
+
         return source;
     }
 

test/freeflow/navierstokes/channel/pipe/main.cc

@@ -101,7 +101,6 @@ int main(int argc, char** argv) try
     SolutionVector sol;
     sol[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
     sol[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
-    problem->applyInitialSolution(sol);
 
     // the grid variables
     using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;

test/freeflow/navierstokes/channel/pipe/params.input

@@ -20,3 +20,6 @@ LiquidDensity = 1000
 
 [Vtk]
 AddVelocity = 1
+
+[Assembly]
+NumericDifference.BaseEpsilon =  1e-3

test/freeflow/navierstokes/channel/pipe/problem.hh

@@ -109,7 +109,7 @@ public:
     { return analyticalSolution(globalPos); }
 
     PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
-    { return analyticalSolution(globalPos); }
+    { return PrimaryVariables(0.0); }
 
     PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
     {
@@ -120,7 +120,7 @@ public:
         const auto y = globalPos[1] - this->gridGeometry().bBoxMin()[1];
         values[Indices::velocityXIdx] = 0.0;
         values[Indices::velocityYIdx] = 2.0*meanInletVelocity_*(1.0 - r*r/(pipeRadius_*pipeRadius_));
-        values[Indices::pressureIdx] = 1e5 + (pipeLength_-y)*meanInletVelocity_*8.0*mu_/(pipeRadius_*pipeRadius_);
+        values[Indices::pressureIdx] = (pipeLength_-y)*meanInletVelocity_*8.0*mu_/(pipeRadius_*pipeRadius_);
 
         return values;
     }