Commit 8572051c authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[freeflow]

added the zeroeq turbulence models from dumux-devel
added zeroeq test problems
to be done:
 - general docu and doxygen docu
 - model description for handbook
 - dune compatibilites 2.3 and >2.3

reviewed by gruenich


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14861 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 09fdce8e
......@@ -40,6 +40,9 @@ AC_CONFIG_FILES([dumux.pc
dumux/freeflow/stokes/Makefile
dumux/freeflow/stokesnc/Makefile
dumux/freeflow/stokesncni/Makefile
dumux/freeflow/zeroeq/Makefile
dumux/freeflow/zeroeqnc/Makefile
dumux/freeflow/zeroeqncni/Makefile
dumux/geomechanics/Makefile
dumux/geomechanics/elastic/Makefile
dumux/geomechanics/el1p2c/Makefile
......@@ -116,6 +119,9 @@ AC_CONFIG_FILES([dumux.pc
test/freeflow/stokes/Makefile
test/freeflow/stokes2c/Makefile
test/freeflow/stokes2cni/Makefile
test/freeflow/zeroeq/Makefile
test/freeflow/zeroeq2c/Makefile
test/freeflow/zeroeq2cni/Makefile
test/geomechanics/Makefile
test/geomechanics/elastic/Makefile
test/geomechanics/el1p2c/Makefile
......
SUBDIRS = stokes stokesnc stokesncni
SUBDIRS = stokes \
stokesnc \
stokesncni \
zeroeq \
zeroeqnc \
zeroeqncni
freeflowdir = $(includedir)/dumux/freeflow
......
zeroeqdir = $(includedir)/dumux/freeflow/zeroeq
zeroeq_HEADERS := $(wildcard *.hh)
include $(top_srcdir)/am/global-rules
// -*- 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
* \brief This file contains the data which is required to calculate
* the fluxes of the ZeroEq model over a face of a finite volume.
*
* This means the methods to calculate the eddy viscosity.
*/
#ifndef DUMUX_ZEROEQ_FLUX_VARIABLES_HH
#define DUMUX_ZEROEQ_FLUX_VARIABLES_HH
#include <dumux/common/math.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/freeflow/stokes/stokesfluxvariables.hh>
#include <dumux/freeflow/zeroeq/zeroeqindices.hh>
#include <dumux/freeflow/zeroeq/zeroeqproperties.hh>
namespace Dumux
{
/*!
* \ingroup ImplicitFluxVariables
* \ingroup BoxZeroEqModel
* \brief This template class contains data which is required to
* calculate the component fluxes over a face of a finite
* volume for a ZeroEq model.
*
* This means the methods to calculate the eddy viscosity.
*/
template <class TypeTag>
class ZeroEqFluxVariables : public GET_PROP_TYPE(TypeTag, BaseStokesFluxVariables)
{
typedef typename GET_PROP_TYPE(TypeTag, BaseStokesFluxVariables) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
enum { dim = GridView::dimension };
typedef typename GridView::template Codim<0>::Entity Element;
typedef Dune::FieldVector<Scalar, dim> DimVector;
public:
ZeroEqFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const int faceIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
: ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
, flowNormal_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, FlowNormal))
, wallNormal_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, WallNormal))
, eddyViscosityModel_(GET_PARAM_FROM_GROUP(TypeTag, int, ZeroEq, EddyViscosityModel))
, karmanConstant_(GET_PROP_VALUE(TypeTag, KarmanConstant))
{
DimVector globalPos = this->face().ipGlobal;
dynamicEddyViscosity_ = 0.0;
mixingLength_ = 0.0;
dynamicEddyViscosityInner_ = 0.0;
dynamicEddyViscosityOuter_ = 0.0;
fz_ = 0.0;
posIdx_ = problem.model().getPosIdx(globalPos);
wallIdx_ = problem.model().getWallIdx(globalPos, posIdx_);
distanceToWallRough_ = std::abs(problem.model().distanceToWallRough(globalPos, wallIdx_, posIdx_));
distanceToWallReal_ = std::abs(problem.model().distanceToWallReal(globalPos, wallIdx_, posIdx_));
for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
maxVelocity_[dimIdx] = problem.model().wall[wallIdx_].maxVelocity[posIdx_][dimIdx];
for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
minVelocity_[dimIdx] = problem.model().wall[wallIdx_].minVelocity[posIdx_][dimIdx];
velGrad_ = this->velocityGrad_[flowNormal_][wallNormal_];
frictionVelocityWall_ = sqrt(problem.model().wall[wallIdx_].wallShearStress[posIdx_]
/ problem.model().wall[wallIdx_].wallDensity[posIdx_]);
yPlusRough_ = distanceToWallRough_ * frictionVelocityWall_ / problem.model().wall[wallIdx_].wallKinematicViscosity[posIdx_];
// calculation of an eddy viscosity only makes sense with Navier-Stokes equation
if (GET_PROP_VALUE(TypeTag, EnableNavierStokes))
calculateEddyViscosity_(problem, element, elemVolVars);
}
protected:
/*!
* \brief This function calculates the dynamic viscosity.
*
* The eddy viscosity is added to the viscosity in stokeslocalresidual.hh at each scv
* face.
*/
void calculateEddyViscosity_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
// no turbulence model
if (eddyViscosityModel_ == EddyViscosityIndices::noEddyViscosityModel)
return;
// Prandtl mixing length
// e.g. Wilcox, D. C., Turbulence Modeling for CFD, 2006
else if (eddyViscosityModel_ == EddyViscosityIndices::prandtl)
{
mixingLength_ = distanceToWallRough_ * karmanConstant_;
dynamicEddyViscosity_ = this->density() * mixingLength() * mixingLength() * std::abs(velGrad_);
}
// modified Van-Driest
// e.g. Bird, Stewart, and Lightfoot, E. N. Transport phenomena, 2007
else if (eddyViscosityModel_ == EddyViscosityIndices::modifiedVanDriest)
{
Scalar aPlus = 26.0;
Scalar bPlus = 0.26;
// eddy viscosity can only be calculated correctly for non-zero distance to walls
mixingLength_ = 0.0;
if (distanceToWallRough_ > 0.0 && yPlusRough_ > 0.0)
mixingLength_= karmanConstant_ * distanceToWallRough_
* (1.0 - std::exp(-yPlusRough_ / aPlus ))
/ std::sqrt(1.0 - std::exp(-bPlus * yPlusRough_));
dynamicEddyViscosity_ = this->density() * mixingLength() * mixingLength() * std::abs(velGrad_);
}
// Baldwin and Lomax
// Baldwin, B. S. & Lomax, H. "Thin Layer Approximation and Algebraic Model for Seperated Turbulent Flows"
// AIAA Journal, 1978, 78--257, 1-9
else if (eddyViscosityModel_ == EddyViscosityIndices::baldwinLomax)
{
// LAW CONSTANTS
const Scalar aPlus = 26.0;
const Scalar cCP = 1.6;
const Scalar cKleb = 0.30;
const Scalar cWK = 0.25;
const Scalar kUpper = 0.0168;
// Calculate muInner
mixingLength_ = 0.0;
if (distanceToWallRough_ > 0.0 && yPlusRough_ > 0.0)
mixingLength_ = karmanConstant_ * distanceToWallRough_
* (1.0 - std::exp(-yPlusRough_ / aPlus ));
Scalar omega1 = this->velocityGrad_[0][1] - this->velocityGrad_[1][0];
Scalar omega2 = 0.0;
Scalar omega3 = 0.0;
if (dim == 3)
{
omega2 = this->velocityGrad_[1][2] - this->velocityGrad_[2][1];
omega3 = this->velocityGrad_[2][0] - this->velocityGrad_[0][2];
}
Scalar omega = sqrt(omega1 * omega1 + omega2 * omega2 + omega3 * omega3);
dynamicEddyViscosityInner_ = this->density() * mixingLength() * mixingLength() * omega;
// Calculate muOuter
fz_ = 0.0;
if (distanceToWallRough_ > 0.0 && yPlusRough_ > 0.0)
fz_ = distanceToWallRough_ * omega * (1.0 - std::exp(-yPlusRough_ / aPlus ));
Scalar fMax = problem.model().wall[wallIdx_].fMax[posIdx_];
Scalar yMax = std::abs(problem.model().wall[wallIdx_].yMax[posIdx_]);
Scalar temp[2] = {0.0, 0.0};
for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
temp[0] += maxVelocity_[dimIdx] * maxVelocity_[dimIdx];
temp[1] += minVelocity_[dimIdx] * minVelocity_[dimIdx];
}
Scalar uDiff = std::sqrt(temp[0]) - std::sqrt(temp[1]);
Scalar f1 = yMax * fMax;
Scalar f2 = cWK * yMax * uDiff * uDiff / fMax;
Scalar fWake = fmin(f1, f2);
Scalar fKleb = 1 / (1 + 5.5 * std::pow(cKleb * distanceToWallRough_ / yMax, 6.0));
dynamicEddyViscosityOuter_ = this->density() * kUpper * cCP * fWake * fKleb;
bool inner = problem.model().useViscosityInner(this->face().ipGlobal, posIdx_);
// muOuter can only be calculated correctly for non-zero fmax
if (fMax == 0)
inner = true;
if (inner)
dynamicEddyViscosity_ = dynamicEddyViscosityInner_;
else
dynamicEddyViscosity_ = dynamicEddyViscosityOuter_;
}
else
{
DUNE_THROW(Dune::NotImplemented, "This eddy viscosity model is not implemented.");
}
Valgrind::CheckDefined(dynamicEddyViscosity_);
Valgrind::CheckDefined(dynamicEddyViscosityInner_);
Valgrind::CheckDefined(dynamicEddyViscosityOuter_);
}
public:
/*!
* \brief Returns the mixing length \f$\mathrm{[m]}\f$ for the element center.
*/
Scalar mixingLength() const
{ return mixingLength_; }
/*!
* \brief Returns the dynamic eddy viscosity
* \f$\mathrm{[Pa \cdot s]} = \mathrm{[N \cdot s/m^2]}\f$.
*/
const Scalar dynamicEddyViscosity() const
{ return dynamicEddyViscosity_; }
/*!
* \brief Returns the kinematic eddy viscosity
* \f$\mathrm{[m^2/s]}\f$ (if implemented).
*/
const Scalar kinematicEddyViscosity() const
{ return dynamicEddyViscosity() / this->density(); }
/*!
* \brief Returns the inner dynamic eddy Viscosity (only Baldwin-Lomax model)
* \f$\mu_\textrm{inner}\f$ in \f$\mathrm{[N\cdot s/m^2]}\f$.
*/
Scalar dynamicEddyViscosityInner() const
{ return dynamicEddyViscosityInner_; }
/*!
* \brief Returns the outer dynamic eddy Viscosity (only Baldwin-Lomax model).
* \f$\mu_\textrm{outer}\f$ in \f$\mathrm{[N\cdot s/m^2]}\f$.
*/
Scalar dynamicEddyViscosityOuter() const
{ return dynamicEddyViscosityOuter_; }
/*!
* \brief Returns the value of the f-function.
*
* \f$\mathrm{f = y \omega \left( 1 - exp \left[ -y^+ / A^+ \right] \right)}\f$.<br>
* y = distanceToWall
*/
Scalar fz() const
{ return fz_; }
/*!
* \brief Returns the friction velocity at the wall \f$\mathrm{[kg/(m*s^2)]}\f$.
*/
Scalar frictionVelocityWall() const
{ return frictionVelocityWall_; }
/*!
* \brief Returns the distance to the corresponding Wall, including surface roughness \f$\mathrm{[m]}\f$.
*/
Scalar distanceToWallRough() const
{ return distanceToWallRough_; }
/*!
* \brief Returns the real distance to the corresponding Wall \f$\mathrm{[m]}\f$.
*/
Scalar distanceToWallReal() const
{ return distanceToWallReal_; }
/*!
* \brief Returns dimensionless wall distance, including surface roughness \f$\mathrm{[-]}\f$.
*/
Scalar yPlusRough() const
{ return yPlusRough_; }
/*!
* \brief Returns the Karman constant \f$\mathrm{[-]}\f$.
*/
Scalar karmanConstant() const
{ return karmanConstant_; }
private:
const int flowNormal_;
const int wallNormal_;
const int eddyViscosityModel_;
const Scalar karmanConstant_;
int wallIdx_;
int posIdx_;
Scalar yPlusRough_;
Scalar distanceToWallReal_;
Scalar distanceToWallRough_;
Scalar sandGrainRoughness_;
Scalar sandGrainRoughnessDimensionless_;
Scalar velGrad_;
Scalar frictionVelocityWall_;
DimVector maxVelocity_;
DimVector minVelocity_;
Scalar mixingLength_;
Scalar dynamicEddyViscosity_;
Scalar dynamicEddyViscosityInner_;
Scalar dynamicEddyViscosityOuter_;
Scalar fz_;
Scalar eps_;
};
} // end namespace
#endif // DUMUX_ZEROEQ_FLUX_VARIABLES_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
* \brief Defines the indices required for the ZeroEq box model.
*/
#ifndef DUMUX_ZEROEQ_INDICES_HH
#define DUMUX_ZEROEQ_INDICES_HH
#include <dumux/freeflow/stokes/stokesindices.hh>
namespace Dumux
{
/*!
* \ingroup BoxZeroEqModel
* \brief The indices for the eddy viscosity model.
*/
struct EddyViscosityIndices
{
// Eddy Viscosity Model Indices
static const int noEddyViscosityModel = 0;
static const int prandtl = 1;
static const int modifiedVanDriest = 2;
static const int baldwinLomax = 3;
};
/*!
* \ingroup BoxZeroEqModel
* \ingroup ImplicitIndices
* \brief The common indices for the isothermal ZeroEq model.
*
* \tparam PVOffset The first index in a primary variable vector.
*/
template <class TypeTag, int PVOffset = 0>
struct ZeroEqCommonIndices : public StokesCommonIndices<TypeTag, PVOffset>
{
static const int scvDataPrecision = 5;
static const int scvDataWidth = scvDataPrecision + 10;
};
} // end namespace
#endif // DUMUX_ZEROEQ_INDICES_HH
This diff is collapsed.
// -*- 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
* \brief Base class for all ZeroEq problems which use the box scheme.
*/
#ifndef DUMUX_ZEROEQ_PROBLEM_HH
#define DUMUX_ZEROEQ_PROBLEM_HH
#include <dumux/freeflow/stokes/stokesmodel.hh>
#include <dumux/freeflow/zeroeq/zeroeqproperties.hh>
namespace Dumux
{
/*!
* \ingroup BoxZeroEqProblems
* \brief Base class for all problems which use the ZeroEq box model.
*
* This implements the call of update functions used for the wall properties of
* the ZeroEq box model.
*/
template<class TypeTag>
class ZeroEqProblem : public StokesProblem<TypeTag>
{
typedef StokesProblem<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
ZeroEqProblem(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
{ }
/*!
* \brief Called by the Dumux::TimeManager in order to initialize the problem.
*
* If you overload this method don't forget to call ParentType::init().<br>
* This initializes all wall-related properties, which are necessary for the
* ZeroEq box model.
*/
void init()
{
// set the initial condition of the model
ParentType::init();
this->model().resetWallProperties();
this->model().resetWallFluidProperties();
}
/*!
* \brief Called by the time manager before the time integration.
*
* This updates all wall-related properties, which are necessary for the
* ZeroEq box model
*/
void preTimeStep()
{
this->model().updateWallProperties();
}
};
}
#endif // DUMUX_ZEROEQ_PROBLEM_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/>. *
*****************************************************************************/
/*!
* \ingroup Properties
* \ingroup ImplicitProperties
* \ingroup BoxZeroEqModel
*
* \file
*
* \brief Defines the supplementary properties required for the
* ZeroEq model.
*
*/
#ifndef DUMUX_ZEROEQ_PROPERTIES_HH
#define DUMUX_ZEROEQ_PROPERTIES_HH
#include <dumux/freeflow/stokes/stokesproperties.hh>
namespace Dumux
{
namespace Properties
{
//////////////////////////////////////////////////////////////////
// Type tags
//////////////////////////////////////////////////////////////////
//! The type tag for the ZeroEq Problem
NEW_TYPE_TAG(BoxZeroEq, INHERITS_FROM(BoxStokes));
//////////////////////////////////////////////////////////////////
// Property tags
//////////////////////////////////////////////////////////////////
NEW_PROP_TAG(KarmanConstant); //!< The Karman constant
NEW_PROP_TAG(BBoxMinIsWall); //!< Sets BBoxMin as a wall
NEW_PROP_TAG(BBoxMaxIsWall); //!< Sets BBoxMax as a wall
NEW_PROP_TAG(ZeroEqFlowNormal); //!< Indicates the main flow direction
NEW_PROP_TAG(ZeroEqWallNormal); //!< Indicates the wall normal direction
NEW_PROP_TAG(ZeroEqBBoxMinSandGrainRoughness); //!< Sets a sand grain roughness at BBoxMin
NEW_PROP_TAG(ZeroEqBBoxMaxSandGrainRoughness); //!< Sets a sand grain roughness at BBoxMax
NEW_PROP_TAG(ZeroEqEddyViscosityModel); //!< Returns used eddy viscosity model
NEW_PROP_TAG(NumberOfIntervals); //!< Returns number of wall intervals
NEW_PROP_TAG(ZeroEqWriteAllSCVData); //!< Returns which scv data are written
NEW_PROP_TAG(BaseStokesModel); //!< Returns the base implementation of the Stokes model
NEW_PROP_TAG(BaseStokesVolumeVariables); //!< Returns the base implementation of the Stokes volume variables
NEW_PROP_TAG(BaseStokesFluxVariables); //!< Returns the base implementation of the Stokes flux variables
}
}
#endif // DUMUX_ZEROEQ_PROPERTIES_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/>. *
*****************************************************************************/
/*!
* \ingroup Properties
* \ingroup ImplicitProperties
* \ingroup BoxZeroEqModel
*
* \file
*
* \brief Defines the properties required for the ZeroEq model.
*/
#ifndef DUMUX_ZEROEQ_PROPERTY_DEFAULTS_HH
#define DUMUX_ZEROEQ_PROPERTY_DEFAULTS_HH
#include "zeroeqfluxvariables.hh"
#include "zeroeqmodel.hh"
namespace Dumux