From f6e09ecf66bed9344fa7506c3bb909abbbd02995 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Mon, 7 Sep 2020 08:06:09 +0200
Subject: [PATCH] [navierstokes][localresidual] Add source term for
 axisymmetric problems

---
 .../navierstokes/staggered/localresidual.hh    | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh
index 64083aa458..6baa841272 100644
--- a/dumux/freeflow/navierstokes/staggered/localresidual.hh
+++ b/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;
     }
 
-- 
GitLab