Commit 220ddc07 authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[rans][lowrekepsilon] Implement first version of a low-Re k-epsilon model

parent 02943424
......@@ -202,6 +202,18 @@
* \brief Zero-equation or algebraic turbulence models
* \copydetails ./freeflow/rans/model.hh
*/
/*!
* \ingroup RANSModel
* \defgroup TwoEqModel 2-Eq. Models
* \brief Two-equation turbulence models
* \copydetails ./freeflow/rans/model.hh
*/
/*!
* \ingroup TwoEqModel
* \defgroup LowReKEpsilonModel Low-Re k-epsilon model
* \brief Low-Re k-epsilon model
* \copydetails ./freeflow/rans/twoeq/lowrekepsilon/model.hh
*/
/*!
* \ingroup FreeflowModels
* \defgroup RANSNCModel Reynolds-Averaged Navier-Stokes nc
......
......@@ -1637,3 +1637,42 @@ url={http://dx.doi.org/10.1007/s11242-015-0599-1}
doi = {10.2514/6.1978-257},
url = {http://dx.doi.org/10.2514/6.1978-257}
}
@Article{Patel1985a,
author = {Patel, V. C. and Rodi, W. and Scheuerer, G.},
title = {{Turbulence models for near-wall and low Reynolds number flows - A review}},
journal = {AIAA Journal},
year = {1985},
volume = {23},
pages = {1308--1319},
number = {9},
doi = {10.2514/3.9086},
url = {http://dx.doi.org/10.2514/3.9086},
}
@Article{Chien1982a,
author = {Chien, Kuei-Yuan},
title = {{Predictions of Channel and Boundary-Layer Flows with a Low-Reynolds-Number Turbulence Model}},
journal = {AIAA Journal},
year = {1982},
volume = {20},
pages = {33--38},
number = {1},
booktitle = {AIAA Journal},
doi = {10.2514/3.51043},
issn = {0001-1452},
publisher = {American Institute of Aeronautics and Astronautics},
url = {http://dx.doi.org/10.2514/3.51043},
}
@Article{Lam1981a,
author = {Lam, C. K. G. and Bremhorst, K. A},
title = {{Modified Form of the k-\varepsilon-Model for Predicting Wall Turbulence}},
journal = {Journal of Fluids Engineering},
year = {1981},
volume = {103},
number = {3},
pages = {456--460},
doi = {10.1115/1.3240815},
url = {http://dx.doi.org/10.1115/1.3240815},
}
add_subdirectory("twoeq")
add_subdirectory("zeroeq")
#install headers
......
......@@ -97,6 +97,7 @@ public:
cellCenter_.resize(this->fvGridGeometry().elementMapper().size(), GlobalPosition(0.0));
velocity_.resize(this->fvGridGeometry().elementMapper().size(), DimVector(0.0));
velocityGradients_.resize(this->fvGridGeometry().elementMapper().size(), DimMatrix(0.0));
stressTensorScalarProduct_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
flowNormalAxis_.resize(this->fvGridGeometry().elementMapper().size(), 0);
wallNormalAxis_.resize(this->fvGridGeometry().elementMapper().size(), 1);
kinematicViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
......@@ -216,7 +217,7 @@ public:
= getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.FlowNormalAxis", -1);
// re-initialize min and max values
velocityMaximum_.assign(this->fvGridGeometry().elementMapper().size(), DimVector(0.0));
velocityMaximum_.assign(this->fvGridGeometry().elementMapper().size(), DimVector(1e-16));
velocityMinimum_.assign(this->fvGridGeometry().elementMapper().size(), DimVector(std::numeric_limits<Scalar>::max()));
// calculate cell-center-averaged velocities
......@@ -283,6 +284,24 @@ public:
{
unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
Dune::FieldMatrix<Scalar, GridView::dimension, GridView::dimension> stressTensor(0.0);
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
{
stressTensor[dimIdx][velIdx] = 0.5 * velocityGradients_[elementID][dimIdx][velIdx]
+ 0.5 * velocityGradients_[elementID][velIdx][dimIdx];
}
}
stressTensorScalarProduct_[elementID] = 0.0;
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
{
stressTensorScalarProduct_[elementID] += stressTensor[dimIdx][velIdx] * stressTensor[dimIdx][velIdx];
}
}
auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
......@@ -338,6 +357,7 @@ public:
std::vector<DimVector> velocityMaximum_;
std::vector<DimVector> velocityMinimum_;
std::vector<DimMatrix> velocityGradients_;
std::vector<Scalar> stressTensorScalarProduct_;
std::vector<unsigned int> wallNormalAxis_;
std::vector<unsigned int> flowNormalAxis_;
std::vector<Scalar> kinematicViscosity_;
......
add_subdirectory("lowrekepsilon")
#install headers
install(FILES
fluxvariables.hh
indices.hh
localresidual.hh
model.hh
# models.hh
problem.hh
volumevariables.hh
vtkoutputfields.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/rans/twoeq/lowrekepsilon)
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup LowReKEpsilonModel
* \copydoc Dumux::LowReKEpsilonFluxVariables
*/
#ifndef DUMUX_LOWREKEPSILON_FLUXVARIABLES_HH
#define DUMUX_LOWREKEPSILON_FLUXVARIABLES_HH
#include <dumux/common/properties.hh>
#include <dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh>
namespace Dumux
{
// forward declaration
template<class TypeTag, DiscretizationMethod discMethod>
class LowReKEpsilonFluxVariablesImpl;
/*!
* \ingroup LowReKEpsilonModel
* \brief The flux variables class for the low-Reynolds k-epsilon model.
This is a convenience alias for that actual,
discretization-specific flux variables.
* \note Not all specializations are currently implemented
*/
template<class TypeTag>
using LowReKEpsilonFluxVariables = LowReKEpsilonFluxVariablesImpl<TypeTag, GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod>;
} // end namespace
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup LowReKEpsilonModel
* \copydoc Dumux::LowReKEpsilonIndices
*/
#ifndef DUMUX_LOWREKEPSILON_INDICES_HH
#define DUMUX_LOWREKEPSILON_INDICES_HH
#include <dumux/freeflow//navierstokes/indices.hh>
namespace Dumux {
// \{
/*!
* \ingroup LowReKEpsilonModel
* \brief The common indices for the isothermal low-Reynolds k-epsilon model.
*
* \tparam dimension The dimension of the problem
*/
template <int dimension>
struct LowReKEpsilonIndices : public NavierStokesIndices<dimension>
{
private:
using ParentType = NavierStokesIndices<dimension>;
public:
static constexpr auto turbulentKineticEnergyEqIdx = ParentType::conti0EqIdx + 1;
static constexpr auto turbulentKineticEnergyIdx = turbulentKineticEnergyEqIdx;
static constexpr auto dissipationEqIdx = turbulentKineticEnergyEqIdx + 1;
static constexpr auto dissipationIdx = dissipationEqIdx;
};
// \}
} // end namespace
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup LowReKEpsilonModel
* \copydoc Dumux::LowReKEpsilonResidual
*/
#ifndef DUMUX_LOWREKEPSILON_LOCAL_RESIDUAL_HH
#define DUMUX_LOWREKEPSILON_LOCAL_RESIDUAL_HH
#include <dumux/common/properties.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/freeflow/navierstokes/localresidual.hh>
#include <dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/localresidual.hh>
namespace Dumux
{
// forward declaration
template<class TypeTag, DiscretizationMethod discMethod>
class LowReKEpsilonResidualImpl;
/*!
* \ingroup LowReKEpsilonModel
* \brief The local residual class for the low-Reynolds k-epsilon model.
This is a convenience alias for the actual,
discretization-specific local residual.
* \note Not all specializations are currently implemented
*/
template<class TypeTag>
using LowReKEpsilonResidual = LowReKEpsilonResidualImpl<TypeTag, GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod>;
}
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup LowReKEpsilonModel
*
* \brief A single-phase, isothermal low-Reynolds k-epsilon model
*
* \copydoc RANSModel
*
* These models calculate the eddy viscosity with two additional PDEs,
* one for the turbulent kinetic energy (k) and for the dissipation (epsilon).
*/
#ifndef DUMUX_LOWREKEPSILON_MODEL_HH
#define DUMUX_LOWREKEPSILON_MODEL_HH
#include <dumux/common/properties.hh>
#include <dumux/freeflow/properties.hh>
#include <dumux/freeflow/rans/model.hh>
#include "fluxvariables.hh"
#include "indices.hh"
#include "localresidual.hh"
#include "volumevariables.hh"
#include "vtkoutputfields.hh"
namespace Dumux
{
namespace Properties {
/*!
* \ingroup LowReKEpsilonModel
* \brief Traits for the low-Reynolds k-epsilon model
*/
template<int dimension>
struct LowReKEpsilonModelTraits
{
//! The dimension of the model
static constexpr int dim() { return dimension; }
//! There are as many momentum balance equations as dimensions,
//! one mass balance equation and two turbulent transport equations
static constexpr int numEq() { return dim()+1+2; }
//! The number of phases is always 1
static constexpr int numPhases() { return 1; }
//! The number of components
static constexpr int numComponents() { return 1; }
//! Enable advection
static constexpr bool enableAdvection() { return true; }
//! The one-phase model has no molecular diffusion
static constexpr bool enableMolecularDiffusion() { return true; }
//! The model is isothermal
static constexpr bool enableEnergyBalance() { return false; }
//! the indices
using Indices = LowReKEpsilonIndices<dim()>;
};
///////////////////////////////////////////////////////////////////////////
// default property values for the isothermal low-Reynolds k-epsilon model
///////////////////////////////////////////////////////////////////////////
//! The type tag for the single-phase, isothermal low-Reynolds k-epsilon model
NEW_TYPE_TAG(LowReKEpsilon, INHERITS_FROM(RANS));
//!< states some specifics of the isothermal low-Reynolds k-epsilon model
SET_PROP(LowReKEpsilon, ModelTraits)
{
private:
using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
static constexpr int dim = GridView::dimension;
public:
using type = LowReKEpsilonModelTraits<dim>;
};
//! The flux variables
SET_TYPE_PROP(LowReKEpsilon, FluxVariables, LowReKEpsilonFluxVariables<TypeTag>);
//! The local residual
SET_TYPE_PROP(LowReKEpsilon, LocalResidual, LowReKEpsilonResidual<TypeTag>);
//! The volume variables
SET_TYPE_PROP(LowReKEpsilon, VolumeVariables, LowReKEpsilonVolumeVariables<TypeTag>);
//! The specific vtk output fields
SET_PROP(LowReKEpsilon, VtkOutputFields)
{
private:
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
public:
using type = LowReKEpsilonVtkOutputFields<FVGridGeometry>;
};
//////////////////////////////////////////////////////////////////
// default property values for the non-isothermal low-Reynolds k-epsilon model
//////////////////////////////////////////////////////////////////
//! The type tag for the single-phase, isothermal low-Reynolds k-epsilon model
NEW_TYPE_TAG(LowReKEpsilonNI, INHERITS_FROM(RANSNI));
//! The volume variables
SET_TYPE_PROP(LowReKEpsilonNI, VolumeVariables, LowReKEpsilonVolumeVariables<TypeTag>);
// \}
}
} // end namespace
#endif // DUMUX_LOWREKEPSILON_MODEL_HH
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup LowReKEpsilonModel
* \copydoc Dumux::LowReKEpsilonModel
*/
#ifndef DUMUX_LOWREKEPSILON_MODELS_HH
#define DUMUX_LOWREKEPSILON_MODELS_HH
namespace Dumux {
/*!
* \ingroup LowReKEpsilonModel
* \brief The available eddy viscosity models
*
* \todo update
* The following models are available:
* -# standard high-Re k-epsilon
* -# Chien, \cite Chien1982a
* -# Lam and Bremhorst, \cite Lam1981a
*
* A good overview and additional models are given in \cite Patel1985a.
*/
class LowReKEpsilonModels
{
public:
static constexpr int standardHighReKEpsilon = 0;
static constexpr int chien = 1;
static constexpr int lamBremhorst = 2;
};
} // end namespace Dumux
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup LowReKEpsilonModel
* \copydoc Dumux::LowReKEpsilonProblem
*/
#ifndef DUMUX_LOWREKEPSILON_PROBLEM_HH
#define DUMUX_LOWREKEPSILON_PROBLEM_HH
#include <dumux/common/properties.hh>
#include <dumux/common/staggeredfvproblem.hh>
#include <dumux/discretization/localview.hh>
#include <dumux/discretization/staggered/elementsolution.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/freeflow/rans/problem.hh>
#include "model.hh"
#include "models.hh"
namespace Dumux
{
/*!
* \ingroup LowReKEpsilonModel
* \brief Low-Re k-epsilon turbulence problem base class.
*
* This implements some base functionality for low-Re k-epsilon models.
*/
template<class TypeTag>
class LowReKEpsilonProblem : public RANSProblem<TypeTag>
{
using ParentType = RANSProblem<TypeTag>;
using Implementation = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
public:
//! The constructor sets the gravity, if desired by the user.
LowReKEpsilonProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
lowReKEpsilonModel_ = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.LowReKEpsilonModel", 1);
if (lowReKEpsilonModel_ != LowReKEpsilonModels::standardHighReKEpsilon
&& lowReKEpsilonModel_ != LowReKEpsilonModels::chien
&& lowReKEpsilonModel_ != LowReKEpsilonModels::lamBremhorst)
{
DUNE_THROW(Dune::NotImplemented, "This low-Re k-epsilon model is not implemented: " << lowReKEpsilonModel_);
}
}
/*!
* \brief Correct size of the static (solution independent) wall variables
*/
void updateStaticWallProperties()
{
ParentType::updateStaticWallProperties();
// update size and initial values of the global vectors
storedDissipationTilde_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
storedKinematicEddyViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
storedTurbulentKineticEnergy_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
}
/*!
* \brief Update the dynamic (solution dependent) relations to the walls
*
* \param curSol The solution vector.
*/
void updateDynamicWallProperties(const SolutionVector& curSol)
{
ParentType::updateDynamicWallProperties(curSol);
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
const int dofIdx = scv.dofIndex();
CellCenterPrimaryVariables priVars(curSol[FVGridGeometry::cellCenterIdx()][dofIdx]);
auto elemSol = elementSolution<FVElementGeometry>(std::move(priVars));
// NOTE: first update the turbulence quantities
storedDissipationTilde_[elementID] = elemSol[0][Indices::dissipationEqIdx];
storedTurbulentKineticEnergy_[elementID] = elemSol[0][Indices::turbulentKineticEnergyEqIdx];
// NOTE: then update the volVars
VolumeVariables volVars;
volVars.update(elemSol, asImp_(), element, scv);
volVars.calculateEddyViscosity();
storedKinematicEddyViscosity_[elementID] = volVars.kinematicEddyViscosity();
}
}
}
/*!
* \brief Returns the index of the used low-Re k-epsilon model