From bc692e2e4d5f11f2d9644dcfd4a2c2bd29051e19 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 30 Mar 2023 19:43:39 +0200
Subject: [PATCH] [exercises][ff-pm][turbulence] Include the SST model and
update the isOnWall stuff
---
.../turbulence/freeflowsubproblem.hh | 102 +++++++++++--
.../models/plotFluxes.py | 2 +-
.../turbulence/freeflowsubproblem.hh | 136 +++++++++++++++---
.../turbulence/properties.hh | 4 +-
4 files changed, 215 insertions(+), 29 deletions(-)
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh b/exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
index 289141a6..5a68f513 100644
--- a/exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
@@ -20,8 +20,8 @@
* \file
* \brief The free-flow sub problem
*/
-#ifndef DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
-#define DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+#ifndef DUMUX_FREEFLOW_TURBULENCE_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW_TURBULENCE_SUBPROBLEM_HH
#include <dumux/common/properties.hh>
#include <dumux/common/boundarytypes.hh>
@@ -46,6 +46,7 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
using ParentType = NavierStokesStaggeredProblem<TypeTag>;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ static constexpr auto dimWorld = GridView::dimensionworld;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
@@ -56,6 +57,7 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
@@ -87,7 +89,23 @@ public:
refTemperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefTemperature");
diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
- getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+ getParamFromGroup<std::string>(this->paramGroup(),
+ "Problem.InterfaceDiffusionCoefficientAvg"));
+
+// TODO: // ******************** uncomment this section for the first task 3.A ****************** //
+// FluidSystem::init();
+// Dumux::TurbulenceProperties<Scalar, dimWorld, true> turbulenceProperties;
+// FluidState fluidState;
+// const auto phaseIdx = 0;
+// fluidState.setPressure(phaseIdx, refPressure_);
+// fluidState.setTemperature(this->spatialParams().temperatureAtPos({}));
+// fluidState.setMassFraction(phaseIdx, phaseIdx, 1.0);
+// Scalar density = FluidSystem::density(fluidState, phaseIdx);
+// Scalar kinematicViscosity = FluidSystem::viscosity(fluidState, phaseIdx) / density;
+// Scalar diameter = this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1];
+// turbulentKineticEnergy_ = turbulenceProperties.turbulentKineticEnergy(refVelocity_, diameter, kinematicViscosity);
+// dissipation_ = turbulenceProperties.dissipationRate(refVelocity_, diameter, kinematicViscosity);
+// // ************************************************************************************* //
}
/*!
@@ -113,7 +131,8 @@ public:
}
// TODO: dumux-course-task 3.A
-// set boundary conditions for the turbulence model at the wall (values.setWall()) at the upper and lower boundary
+// set wall conditions for the turbulence model at the wall (values.setWall()) at the upper and lower boundary
+// set boundary conditions for the turbulence model primary variables everywhere (outflow on right boundary, otherwise dirichlet)
if (onLowerBoundary_(globalPos))
{
values.setDirichlet(Indices::velocityXIdx);
@@ -153,16 +172,65 @@ public:
}
/*!
- * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ * \brief Evaluate the boundary conditions for dirichlet values at the boundary.
*
- * \param element The element
- * \param scvf The subcontrolvolume face
+ * \param element The finite element
+ * \param scvf the sub control volume face
+ * \note used for cell-centered discretization schemes
*/
- PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+ PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
{
- PrimaryVariables values(0.0);
- values = initialAtPos(pos);
+ const auto globalPos = scvf.ipGlobal();
+ PrimaryVariables values(initialAtPos(globalPos));
+ return values;
+ }
+
+ /*!
+ * \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param scv The sub control volume
+ * \param pvIdx The primary variable index in the solution vector
+ */
+ template<class Element, class FVElementGeometry, class SubControlVolume>
+ bool isDirichletCell(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolume& scv,
+ int pvIdx) const
+ {
+ if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::sst)
+ {
+ // For the komega model we set a fixed dissipation (omega) for all cells at the wall
+ for (const auto& scvf : scvfs(fvGeometry))
+ if (this->boundaryTypes(element, scvf).hasWall() && pvIdx == Indices::dissipationIdx)
+ return true;
+ }
+ return false;
+ }
+ /*!
+ * \brief Evaluate the boundary conditions for fixed values at cell centers
+ *
+ * \param element The finite element
+ * \param scv the sub control volume
+ * \note used for cell-centered discretization schemes
+ */
+ PrimaryVariables dirichlet([[maybe_unused]] const Element& element, const SubControlVolume& scv) const
+ {
+ const auto globalPos = scv.center();
+ PrimaryVariables values(initialAtPos(globalPos));
+ unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
+
+ if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::sst)
+ {
+ // For the komega model we set a fixed value for the dissipation
+ const auto wallDistance = ParentType::wallDistance(elementIdx);
+ using std::pow;
+ values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity(elementIdx)
+ / (ParentType::betaOmega() * wallDistance * wallDistance);
+ return values;
+ }
return values;
}
@@ -220,6 +288,18 @@ public:
values[Indices::velocityXIdx] = refVelocity();
values[Indices::temperatureIdx] = refTemperature();
+ // TODO: dumux-course-task 3.A
+ // Set initial conditions for the TKE and the Dissipation. Values calculated in the constructor
+// values[Indices::turbulentKineticEnergyIdx] = TODO??;
+// values[Indices::dissipationIdx] = TODO??;
+// // TODO: dumux-course-task 3.B
+// // Remove the condition `onUpperBoundary_(globalPos)` here.
+// if(onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
+// {
+// values[Indices::turbulentKineticEnergyIdx] = 0.0;
+// values[Indices::dissipationIdx] = 0.0;
+// }
+
// TODO: dumux-course-task 3.B
// Remove the condition `onUpperBoundary_(globalPos)` here.
if(onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
@@ -317,6 +397,8 @@ private:
Scalar refPressure_;
Scalar refMoleFrac_;
Scalar refTemperature_;
+ Scalar turbulentKineticEnergy_;
+ Scalar dissipation_;
TimeLoopPtr timeLoop_;
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/plotFluxes.py b/exercises/solution/exercise-coupling-ff-pm/models/plotFluxes.py
index a3955cee..667b97bf 120000
--- a/exercises/solution/exercise-coupling-ff-pm/models/plotFluxes.py
+++ b/exercises/solution/exercise-coupling-ff-pm/models/plotFluxes.py
@@ -1 +1 @@
-/home/nedc/DUMUX_summerschool/dumux-course/exercises/exercise-coupling-ff-pm/models/plotFluxes.py
\ No newline at end of file
+../../../exercise-coupling-ff-pm/models/plotFluxes.py
\ No newline at end of file
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh b/exercises/solution/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
index 41fae5c3..66dbbc3b 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
@@ -20,22 +20,25 @@
* \file
* \brief The free-flow sub problem
*/
-#ifndef DUMUX_FREEFLOW_TURBULENCE_SUBPROBLEM_HH
-#define DUMUX_FREEFLOW_TURBULENCE_SUBPROBLEM_HH
+#ifndef DUMUX_FREEFLOW_TURBULENCE_SUBPROBLEM_SOL_HH
+#define DUMUX_FREEFLOW_TURBULENCE_SUBPROBLEM_SOL_HH
-#if EXNUMBER >= 1
-#include <dumux/freeflow/rans/problem.hh>
-#include <dumux/freeflow/rans/boundarytypes.hh>
-#else
-#include <dumux/freeflow/navierstokes/staggered/problem.hh>
-#endif
#include <dumux/common/timeloop.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/multidomain/boundary/stokesdarcy/couplingdata.hh>
+
+#if EXNUMBER >= 1
+#include <dumux/freeflow/turbulencemodel.hh>
+#include <dumux/freeflow/turbulenceproperties.hh>
+#include <dumux/freeflow/rans/twoeq/sst/problem.hh>
+#include <dumux/freeflow/rans/boundarytypes.hh>
+#else
+#include <dumux/freeflow/navierstokes/staggered/problem.hh>
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
+#endif
namespace Dumux {
@@ -54,6 +57,7 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
#endif
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+ static constexpr auto dimWorld = GridView::dimensionworld;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
@@ -61,11 +65,11 @@ class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
#if EXNUMBER >= 1
using BoundaryTypes = Dumux::RANSBoundaryTypes<ModelTraits, ModelTraits::numEq()>;
#else
- using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag,
- Properties::ModelTraits>::numEq()>;
+ using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
#endif
using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
@@ -98,6 +102,21 @@ public:
diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+#if EXNUMBER >=1
+ FluidSystem::init();
+ Dumux::TurbulenceProperties<Scalar, dimWorld, true> turbulenceProperties;
+ FluidState fluidState;
+ const auto phaseIdx = 0;
+ fluidState.setPressure(phaseIdx, refPressure_);
+ fluidState.setTemperature(this->spatialParams().temperatureAtPos({}));
+ fluidState.setMassFraction(phaseIdx, phaseIdx, refMoleFrac_);
+ Scalar density = FluidSystem::density(fluidState, phaseIdx);
+ Scalar kinematicViscosity = FluidSystem::viscosity(fluidState, phaseIdx) / density;
+ Scalar diameter = this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1];
+ turbulentKineticEnergy_ = turbulenceProperties.turbulentKineticEnergy(refVelocity_, diameter, kinematicViscosity);
+ dissipation_ = turbulenceProperties.dissipationRate(refVelocity_, diameter, kinematicViscosity);
+#endif
+
}
/*!
@@ -114,6 +133,20 @@ public:
const auto& globalPos = scvf.center();
+#if EXNUMBER >=1
+ if(onRightBoundary_(globalPos))
+ {
+ values.setOutflow(Indices::turbulentKineticEnergyEqIdx);
+ values.setOutflow(Indices::dissipationEqIdx);
+ }
+ else
+ {
+ // walls and inflow
+ values.setDirichlet(Indices::turbulentKineticEnergyIdx);
+ values.setDirichlet(Indices::dissipationIdx);
+ }
+#endif
+
if (onLeftBoundary_(globalPos))
{
values.setDirichlet(Indices::velocityXIdx);
@@ -169,20 +202,70 @@ public:
values.setCouplingNeumann(Indices::momentumYBalanceIdx);
values.setBeaversJoseph(Indices::momentumXBalanceIdx);
}
+
return values;
}
/*!
- * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ * \brief Evaluate the boundary conditions for a dirichlet values at the boundary.
*
- * \param element The element
- * \param scvf The subcontrolvolume face
+ * \param element The finite element
+ * \param scvf the sub control volume face
+ * \note used for cell-centered discretization schemes
*/
- PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+ PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
{
- PrimaryVariables values(0.0);
- values = initialAtPos(pos);
+ const auto globalPos = scvf.ipGlobal();
+ PrimaryVariables values(initialAtPos(globalPos));
+ return values;
+ }
+ /*!
+ * \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param scv The sub control volume
+ * \param pvIdx The primary variable index in the solution vector
+ */
+ template<class Element, class FVElementGeometry, class SubControlVolume>
+ bool isDirichletCell(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolume& scv,
+ int pvIdx) const
+ {
+ if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::sst)
+ {
+ // For the komega model we set a fixed dissipation (omega) for all cells at the wall
+ for (const auto& scvf : scvfs(fvGeometry))
+ if (this->boundaryTypes(element, scvf).hasWall() && pvIdx == Indices::dissipationIdx)
+ return true;
+ }
+ return false;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for fixed values at cell centers
+ *
+ * \param element The finite element
+ * \param scv the sub control volume
+ * \note used for cell-centered discretization schemes
+ */
+ PrimaryVariables dirichlet([[maybe_unused]] const Element& element, const SubControlVolume& scv) const
+ {
+ const auto globalPos = scv.center();
+ PrimaryVariables values(initialAtPos(globalPos));
+ unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
+
+ if constexpr (ModelTraits::turbulenceModel() == TurbulenceModel::sst)
+ {
+ // For the komega model we set a fixed value for the dissipation
+ const auto wallDistance = ParentType::wallDistance(elementIdx);
+ using std::pow;
+ values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity(elementIdx)
+ / (ParentType::betaOmega() * wallDistance * wallDistance);
+ return values;
+ }
return values;
}
@@ -240,6 +323,25 @@ public:
values[Indices::velocityXIdx] = refVelocity();
values[Indices::temperatureIdx] = refTemperature();
+#if EXNUMBER >= 1
+ values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_;
+ values[Indices::dissipationIdx] = dissipation_;
+#endif
+
+#if EXNUMBER == 1
+ if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos))
+ {
+ values[Indices::turbulentKineticEnergyIdx] = 0.0;
+ values[Indices::dissipationIdx] = 0.0;
+ }
+#elif EXNUMBER >= 2
+ if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos))
+ {
+ values[Indices::turbulentKineticEnergyIdx] = 0.0;
+ values[Indices::dissipationIdx] = 0.0;
+ }
+#endif
+
#if EXNUMBER >= 2
if(onLowerBoundary_(globalPos))
values[Indices::velocityXIdx] = 0.0;
@@ -340,6 +442,8 @@ private:
Scalar refPressure_;
Scalar refMoleFrac_;
Scalar refTemperature_;
+ Scalar turbulentKineticEnergy_;
+ Scalar dissipation_;
TimeLoopPtr timeLoop_;
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/properties.hh b/exercises/solution/exercise-coupling-ff-pm/turbulence/properties.hh
index 5389371e..3b6b962a 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/properties.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/properties.hh
@@ -43,7 +43,7 @@
#include <dumux/discretization/staggered/freeflow/properties.hh>
#if EXNUMBER >= 1
-#include <dumux/freeflow/compositional/zeroeqncmodel.hh>
+#include <dumux/freeflow/compositional/sstncmodel.hh>
#else
#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
#endif
@@ -56,7 +56,7 @@ namespace Dumux::Properties {
namespace TTag {
struct PorousMediumModel { using InheritsFrom = std::tuple<TwoPTwoCNI, CCTpfaModel>; };
#if EXNUMBER >= 1
-struct FreeflowModel { using InheritsFrom = std::tuple<ZeroEqNCNI, StaggeredFreeFlowModel>; };
+struct FreeflowModel { using InheritsFrom = std::tuple<SSTNCNI, StaggeredFreeFlowModel>; };
#else
struct FreeflowModel { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
#endif
--
GitLab