Skip to content
Snippets Groups Projects
Commit 80bfcd7f authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[freeflow][komega] Clean-up

parent c11cee90
No related branches found
No related tags found
1 merge request!1027Freeflow/komega
...@@ -104,8 +104,6 @@ public: ...@@ -104,8 +104,6 @@ public:
= ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermK); = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermK);
flux[dissipationEqIdx] flux[dissipationEqIdx]
= ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermOmega); = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermOmega);
Dune::dverb << " k_adv " << ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermK);
Dune::dverb << " w_adv " << ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermOmega);
// calculate diffusive flux // calculate diffusive flux
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
...@@ -154,9 +152,6 @@ public: ...@@ -154,9 +152,6 @@ public:
+= coeff_k / distance += coeff_k / distance
* (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy()) * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy())
* scvf.area(); * scvf.area();
Dune::dverb << " k_diff " << coeff_k / distance
* (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy())
* scvf.area();
} }
if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::dissipationEqIdx) if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::dissipationEqIdx)
|| bcTypes.isSymmetry()))) || bcTypes.isSymmetry())))
...@@ -165,9 +160,6 @@ public: ...@@ -165,9 +160,6 @@ public:
+= coeff_w / distance += coeff_w / distance
* (insideVolVars.dissipation() - outsideVolVars.dissipation()) * (insideVolVars.dissipation() - outsideVolVars.dissipation())
* scvf.area(); * scvf.area();
Dune::dverb << " w_diff " << coeff_w / distance
* (insideVolVars.dissipation() - outsideVolVars.dissipation())
* scvf.area();
} }
return flux; return flux;
} }
......
...@@ -99,21 +99,22 @@ public: ...@@ -99,21 +99,22 @@ public:
CellCenterPrimaryVariables source = ParentType::computeSourceForCellCenter(problem, element, fvGeometry, CellCenterPrimaryVariables source = ParentType::computeSourceForCellCenter(problem, element, fvGeometry,
elemVolVars, elemFaceVars, scv); elemVolVars, elemFaceVars, scv);
using std::min;
const auto& volVars = elemVolVars[scv]; const auto& volVars = elemVolVars[scv];
// production // production
Scalar productionTerm = 2.0 * volVars.kinematicEddyViscosity()* volVars.stressTensorScalarProduct(); Scalar productionTerm = 2.0 * volVars.kinematicEddyViscosity() * volVars.stressTensorScalarProduct();
if (ModelTraits::enableKOmegaProductionLimiter()) if (ModelTraits::enableKOmegaProductionLimiter())
{ {
Scalar productionAlternative = 20.0 * volVars.betaK() * volVars.turbulentKineticEnergy() * volVars.dissipation(); Scalar productionAlternative = 20.0 * volVars.betaK() * volVars.turbulentKineticEnergy() * volVars.dissipation();
productionTerm = std::min(productionTerm, productionAlternative); productionTerm = min(productionTerm, productionAlternative);
} }
source[turbulentKineticEnergyEqIdx] += productionTerm; source[turbulentKineticEnergyEqIdx] += productionTerm;
source[dissipationEqIdx] += volVars.alpha() * ( volVars.dissipation() / volVars.turbulentKineticEnergy() ) * productionTerm; source[dissipationEqIdx] += volVars.alpha() * volVars.dissipation() / volVars.turbulentKineticEnergy() * productionTerm;
// destruction // destruction
source[turbulentKineticEnergyEqIdx] -= volVars.betaK() * volVars.turbulentKineticEnergy() * volVars.dissipation(); source[turbulentKineticEnergyEqIdx] -= volVars.betaK() * volVars.turbulentKineticEnergy() * volVars.dissipation();
source[dissipationEqIdx] -= volVars.betaOmega() * volVars.dissipation() * volVars.dissipation(); source[dissipationEqIdx] -= volVars.betaOmega() * volVars.dissipation() * volVars.dissipation();
return source; return source;
} }
......
...@@ -98,7 +98,6 @@ public: ...@@ -98,7 +98,6 @@ public:
storedDissipation_ = problem.storedDissipation_[RANSParentType::elementID()]; storedDissipation_ = problem.storedDissipation_[RANSParentType::elementID()];
storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementID()]; storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementID()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementID()]; stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementID()];
dofPosition_ = scv.dofPosition();
if (problem.useStoredEddyViscosity_) if (problem.useStoredEddyViscosity_)
dynamicEddyViscosity_ = problem.storedDynamicEddyViscosity_[RANSParentType::elementID()]; dynamicEddyViscosity_ = problem.storedDynamicEddyViscosity_[RANSParentType::elementID()];
else else
...@@ -200,15 +199,15 @@ public: ...@@ -200,15 +199,15 @@ public:
//! \brief Returns the \$f \alpha \$f value //! \brief Returns the \$f \alpha \$f value
const Scalar alpha() const const Scalar alpha() const
{ return 0.520; } { return 0.52; }
//! \brief Returns the \$f \sigma_k \$f constant //! \brief Returns the \$f \sigma_k \$f constant
const Scalar sigmaK() const const Scalar sigmaK() const
{ return 0.60; } { return 0.6; }
//! \brief Returns the \$f \sigma_{\omega} \$f constant //! \brief Returns the \$f \sigma_{\omega} \$f constant
const Scalar sigmaOmega() const const Scalar sigmaOmega() const
{ return 0.50; } { return 0.5; }
//! \brief Returns the \$f \beta_k \$f constant //! \brief Returns the \$f \beta_k \$f constant
const Scalar betaK() const const Scalar betaK() const
...@@ -237,7 +236,6 @@ public: ...@@ -237,7 +236,6 @@ public:
protected: protected:
Scalar betaOmega_; Scalar betaOmega_;
Dune::FieldVector<Scalar, Traits::ModelTraits::dim()> dofPosition_;
Scalar eddyDiffusivity_; Scalar eddyDiffusivity_;
Scalar dynamicEddyViscosity_; Scalar dynamicEddyViscosity_;
Scalar dissipation_; Scalar dissipation_;
......
...@@ -51,10 +51,8 @@ public: ...@@ -51,10 +51,8 @@ public:
template <class VtkOutputModule> template <class VtkOutputModule>
static void add(VtkOutputModule& vtk) static void add(VtkOutputModule& vtk)
{ {
vtk.addVolumeVariable([](const auto& v){ return v.turbulentKineticEnergy(); }, "turbulentKineticEnergy"); vtk.addVolumeVariable([](const auto& v){ return v.turbulentKineticEnergy(); }, "k");
vtk.addVolumeVariable([](const auto& v){ return v.dissipation(); }, "dissipation"); vtk.addVolumeVariable([](const auto& v){ return v.dissipation(); }, "omega");
vtk.addVolumeVariable([](const auto& v){ return v.stressTensorScalarProduct(); }, "stressTensorScalarProduct");
vtk.addVolumeVariable([](const auto& v){ return v.kinematicEddyViscosity(); }, "eddyViscosity");
vtk.addVolumeVariable([](const auto& v){ return 2.0 * v.kinematicEddyViscosity() * v.stressTensorScalarProduct();}, "production_k"); vtk.addVolumeVariable([](const auto& v){ return 2.0 * v.kinematicEddyViscosity() * v.stressTensorScalarProduct();}, "production_k");
vtk.addVolumeVariable([](const auto& v){ return 2.0 * v.kinematicEddyViscosity() * v.stressTensorScalarProduct() * v.alpha() * ( v.dissipation() / v.turbulentKineticEnergy() ) ; }, "production_omega"); vtk.addVolumeVariable([](const auto& v){ return 2.0 * v.kinematicEddyViscosity() * v.stressTensorScalarProduct() * v.alpha() * ( v.dissipation() / v.turbulentKineticEnergy() ) ; }, "production_omega");
vtk.addVolumeVariable([](const auto& v){ return v.betaK() * v.turbulentKineticEnergy() * v.dissipation() ;}, "destruction_k"); vtk.addVolumeVariable([](const auto& v){ return v.betaK() * v.turbulentKineticEnergy() * v.dissipation() ;}, "destruction_k");
......
...@@ -274,7 +274,7 @@ int main(int argc, char** argv) try ...@@ -274,7 +274,7 @@ int main(int argc, char** argv) try
gnuplot_velocityProfile.addFileToPlot(std::string(fileName) + ".csv", "u 7:($28/0.2456) w l lc 7"); gnuplot_velocityProfile.addFileToPlot(std::string(fileName) + ".csv", "u 7:($28/0.2456) w l lc 7");
#elif KOMEGA #elif KOMEGA
gnuplot_velocityProfile.addFileToPlot("pdelab-komega.csv", "u 5:($29/0.2456) w l lw 2 t 'PDELab k-omega'"); gnuplot_velocityProfile.addFileToPlot("pdelab-komega.csv", "u 5:($29/0.2456) w l lw 2 t 'PDELab k-omega'");
gnuplot_velocityProfile.addFileToPlot(std::string(fileName) + ".csv", "u 7:($32/0.2456) w l lc 7"); gnuplot_velocityProfile.addFileToPlot(std::string(fileName) + ".csv", "u 7:($30/0.2456) w l lc 7");
#else #else
gnuplot_velocityProfile.addFileToPlot("pdelab-zeroeq.csv", "u 5:($26/0.2456) w l lw 2 t 'PDELab 0-Eq.'"); gnuplot_velocityProfile.addFileToPlot("pdelab-zeroeq.csv", "u 5:($26/0.2456) w l lw 2 t 'PDELab 0-Eq.'");
gnuplot_velocityProfile.addFileToPlot(std::string(fileName) + ".csv", "u 7:($24/0.2456) w l lc 7"); gnuplot_velocityProfile.addFileToPlot(std::string(fileName) + ".csv", "u 7:($24/0.2456) w l lc 7");
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment