From c11cee90bc3bac6f8037433382a155e0052687ce Mon Sep 17 00:00:00 2001
From: Thomas Fetzer <thomas.fetzer@iws.uni-stuttgart.de>
Date: Fri, 22 Jun 2018 08:54:28 +0200
Subject: [PATCH] [freeflow][komega] Implement fixed-value cells at walls and
improve BCs
---
.../twoeq/komega/staggered/localresidual.hh | 20 +++++++++
test/freeflow/rans/pipelauferproblem.hh | 44 ++++++++-----------
2 files changed, 38 insertions(+), 26 deletions(-)
diff --git a/dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh b/dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh
index 034e9754d8..baaa228d55 100644
--- a/dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh
+++ b/dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh
@@ -65,6 +65,7 @@ class KOmegaResidualImpl<TypeTag, BaseLocalResidual, DiscretizationMethod::stagg
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using CellCenterResidual = CellCenterPrimaryVariables;
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
@@ -116,6 +117,25 @@ public:
return source;
}
+
+ /*!
+ * \brief Sets a fixed Dirichlet value for a cell (such as for dissipation) at the boundary.
+ */
+ template<class BoundaryTypes>
+ void setFixedCell(CellCenterResidual& residual,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& insideScv,
+ const ElementVolumeVariables& elemVolVars,
+ const BoundaryTypes& bcTypes) const
+ {
+ // set a fixed dissipation for cells adjacent to a wall
+ if(bcTypes.isDirichletCell(Indices::dissipationEqIdx))
+ {
+ const auto& insideVolVars = elemVolVars[insideScv];
+ residual[dissipationEqIdx] = insideVolVars.dissipation() - problem.dirichlet(element, insideScv)[Indices::dissipationEqIdx];
+ }
+ }
};
}
diff --git a/test/freeflow/rans/pipelauferproblem.hh b/test/freeflow/rans/pipelauferproblem.hh
index e2c9f2ead6..2442b49f04 100644
--- a/test/freeflow/rans/pipelauferproblem.hh
+++ b/test/freeflow/rans/pipelauferproblem.hh
@@ -153,8 +153,11 @@ public:
Scalar diameter = this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1];
viscosityTilde_ = turbulenceProperties.viscosityTilde(inletVelocity_, diameter, kinematicViscosity);
turbulentKineticEnergy_ = turbulenceProperties.turbulentKineticEnergy(inletVelocity_, diameter, kinematicViscosity);
+#if KOMEGA
+ dissipation_ = turbulenceProperties.dissipationRate(inletVelocity_, diameter, kinematicViscosity);
+#else
dissipation_ = turbulenceProperties.dissipation(inletVelocity_, diameter, kinematicViscosity);
- dissipationRate_ = turbulenceProperties.dissipationRate(inletVelocity_, diameter, kinematicViscosity);
+#endif
std::cout << std::endl;
}
@@ -235,14 +238,19 @@ public:
#if LOWREKEPSILON || KEPSILON || KOMEGA
values.setDirichlet(Indices::turbulentKineticEnergyIdx);
values.setDirichlet(Indices::dissipationIdx);
+
+#if KOMEGA
+ // set a fixed dissipation (omega) in one cell
+ if (isOnWall(globalPos))
+ values.setDirichletCell(Indices::dissipationIdx);
+#endif
#endif
}
return values;
}
/*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * control volume face.
+ * \brief Evaluate the boundary conditions for a dirichlet values at the boundary.
*
* \param element The finite element
* \param scvf the sub control volume face
@@ -252,11 +260,6 @@ public:
{
const auto globalPos = scvf.ipGlobal();
PrimaryVariables values(initialAtPos(globalPos));
-#if KOMEGA
- unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
- const auto wallDistance = ParentType::wallDistance_[elementID];
- values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementID] / (ParentType::betaOmega() * std::pow(wallDistance, 2));
-#endif
#if NONISOTHERMAL
if (time() > 10.0)
{
@@ -271,8 +274,7 @@ public:
}
/*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * control volume center.
+ * \brief Evaluate the boundary conditions for fixed values at cell centers
*
* \param element The finite element
* \param scv the sub control volume
@@ -283,9 +285,11 @@ public:
const auto globalPos = scv.center();
PrimaryVariables values(initialAtPos(globalPos));
#if KOMEGA
-// unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
-// const auto wallDistance = ParentType::wallDistance_[elementID];
-// values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementID] / (ParentType::betaOmega() * std::pow(wallDistance, 2));
+ using std::pow;
+ unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+ const auto wallDistance = ParentType::wallDistance_[elementID];
+ values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementID]
+ / (ParentType::betaOmega() * pow(wallDistance, 2));
#endif
return values;
}
@@ -314,7 +318,7 @@ public:
}
#endif
-#if LOWREKEPSILON || KEPSILON
+#if LOWREKEPSILON || KEPSILON || KOMEGA
values[Indices::turbulentKineticEnergyEqIdx] = turbulentKineticEnergy_;
values[Indices::dissipationEqIdx] = dissipation_;
if (isOnWall(globalPos))
@@ -322,17 +326,6 @@ public:
values[Indices::turbulentKineticEnergyEqIdx] = 0.0;
values[Indices::dissipationEqIdx] = 0.0;
}
-#elif KOMEGA
- if (time() < eps_ && startWithZeroVelocity_)
- {
- values[Indices::velocityXIdx] = 0.0;
- }
- values[Indices::turbulentKineticEnergyEqIdx] = turbulentKineticEnergy_;
- values[Indices::dissipationEqIdx] = dissipation_;
- if (isOnWall(globalPos))
- {
- values[Indices::turbulentKineticEnergyEqIdx] = 0.0;
- }
#endif
return values;
@@ -370,7 +363,6 @@ private:
Scalar viscosityTilde_;
Scalar turbulentKineticEnergy_;
Scalar dissipation_;
- Scalar dissipationRate_;
TimeLoopPtr timeLoop_;
};
} //end namespace
--
GitLab