Commit 63fd0cb2 authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[lowrekepsilon] Adapat to changes in 985e5cf2

parent c862561e
......@@ -101,8 +101,20 @@ 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>);
//! Set the volume variables property
SET_PROP(LowReKEpsilon, VolumeVariables)
{
private:
using PV = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FSY = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using FST = typename GET_PROP_TYPE(TypeTag, FluidState);
using MT = typename GET_PROP_TYPE(TypeTag, ModelTraits);
using Traits = NavierStokesVolumeVariablesTraits<PV, FSY, FST, MT>;
using NSVolVars = NavierStokesVolumeVariables<Traits>;
public:
using type = LowReKEpsilonVolumeVariables<Traits, NSVolVars>;
};
//! The specific vtk output fields
SET_PROP(LowReKEpsilon, VtkOutputFields)
......@@ -120,8 +132,20 @@ public:
//! 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>);
//! Set the volume variables property
SET_PROP(LowReKEpsilonNI, VolumeVariables)
{
private:
using PV = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FSY = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using FST = typename GET_PROP_TYPE(TypeTag, FluidState);
using MT = typename GET_PROP_TYPE(TypeTag, ModelTraits);
using Traits = NavierStokesVolumeVariablesTraits<PV, FSY, FST, MT>;
using NSVolVars = NavierStokesVolumeVariables<Traits>;
public:
using type = LowReKEpsilonVolumeVariables<Traits, NSVolVars>;
};
// \}
}
......
......@@ -85,7 +85,7 @@ public:
// 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);
storedDynamicEddyViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
storedTurbulentKineticEnergy_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
}
......@@ -116,8 +116,7 @@ public:
// NOTE: then update the volVars
VolumeVariables volVars;
volVars.update(elemSol, asImp_(), element, scv);
volVars.calculateEddyViscosity();
storedKinematicEddyViscosity_[elementID] = volVars.kinematicEddyViscosity();
storedDynamicEddyViscosity_[elementID] = volVars.calculateEddyViscosity();
}
}
}
......@@ -130,7 +129,7 @@ public:
public:
std::vector<Scalar> storedDissipationTilde_;
std::vector<Scalar> storedKinematicEddyViscosity_;
std::vector<Scalar> storedDynamicEddyViscosity_;
std::vector<Scalar> storedTurbulentKineticEnergy_;
int lowReKEpsilonModel_;
bool useStoredEddyViscosity_;
......
......@@ -35,38 +35,27 @@
namespace Dumux
{
// forward declaration
template <class TypeTag, bool enableEnergyBalance>
class LowReKEpsilonVolumeVariablesImplementation;
/*!
* \ingroup LowReKEpsilonModel
* \brief Volume variables for the single-phase 0-Eq. model.
* The class is specialized for isothermal and non-isothermal models.
*/
template <class TypeTag>
using LowReKEpsilonVolumeVariables = LowReKEpsilonVolumeVariablesImplementation<TypeTag, GET_PROP_TYPE(TypeTag, ModelTraits)::enableEnergyBalance()>;
/*!
* \ingroup LowReKEpsilonModel
* \brief Volume variables for the isothermal single-phase 0-Eq. model.
*/
template <class TypeTag>
class LowReKEpsilonVolumeVariablesImplementation<TypeTag, false>
: virtual public RANSVolumeVariablesImplementation<TypeTag, false>
template <class Traits, class NSVolumeVariables>
class LowReKEpsilonVolumeVariables
: public RANSVolumeVariables< Traits, LowReKEpsilonVolumeVariables<Traits, NSVolumeVariables> >
, public NSVolumeVariables
{
using ParentType = RANSVolumeVariablesImplementation<TypeTag, false>;
using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using ThisType = LowReKEpsilonVolumeVariables<Traits, NSVolumeVariables>;
using RANSParentType = RANSVolumeVariables<Traits, ThisType>;
using NavierStokesParentType = NSVolumeVariables;
using Scalar = typename Traits::PrimaryVariables::value_type;
using Indices = typename Traits::ModelTraits::Indices;
static constexpr bool enableEnergyBalance = Traits::ModelTraits::enableEnergyBalance();
public:
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
//! export the underlying fluid system
using FluidSystem = typename Traits::FluidSystem;
/*!
* \brief Update all quantities for a given control volume
......@@ -77,13 +66,13 @@ public:
* \param element An element which contains part of the control volume
* \param scv The sub-control volume
*/
template<class ElementSolution>
template<class ElementSolution, class Problem, class Element, class SubControlVolume>
void update(const ElementSolution &elemSol,
const Problem &problem,
const Element &element,
const SubControlVolume& scv)
{
ParentType::update(elemSol, problem, element, scv);
NavierStokesParentType::update(elemSol, problem, element, scv);
updateRANSProperties(elemSol, problem, element, scv);
}
......@@ -97,44 +86,56 @@ public:
* \param element An element which contains part of the control volume
* \param scv The sub-control volume
*/
template<class ElementSolution>
template<class ElementSolution, class Problem, class Element, class SubControlVolume>
void updateRANSProperties(const ElementSolution &elemSol,
const Problem &problem,
const Element &element,
const SubControlVolume& scv)
{
ParentType::updateRANSProperties(elemSol, problem, element, scv);
RANSParentType::updateRANSProperties(elemSol, problem, element, scv);
lowReKEpsilonModel_ = problem.lowReKEpsilonModel();
turbulentKineticEnergy_ = elemSol[0][Indices::turbulentKineticEnergyIdx];
dissipationTilde_ = elemSol[0][Indices::dissipationIdx];
storedDissipationTilde_ = problem.storedDissipationTilde_[ParentType::elementID()];
storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[ParentType::elementID()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[ParentType::elementID()];
storedDissipationTilde_ = problem.storedDissipationTilde_[RANSParentType::elementID()];
storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementID()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementID()];
if (problem.useStoredEddyViscosity_)
ParentType::setKinematicEddyViscosity(problem.storedKinematicEddyViscosity_[ParentType::elementID()]);
dynamicEddyViscosity_ = problem.storedDynamicEddyViscosity_[RANSParentType::elementID()];
else
calculateEddyViscosity();
dynamicEddyViscosity_ = calculateEddyViscosity();
}
/*!
* \brief Return the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$ of the flow
*/
Scalar dynamicEddyViscosity() const
{ return dynamicEddyViscosity_; }
/*!
* \brief Return the effective dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
* control volume.
*/
Scalar effectiveViscosity() const
{ return NavierStokesParentType::viscosity() + dynamicEddyViscosity(); }
/*!
* \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
* of the fluid-flow in the sub-control volume.
*/
template<bool eB = enableEnergyBalance, typename std::enable_if_t<eB, int> = 0>
Scalar effectiveThermalConductivity() const
{
return NavierStokesParentType::thermalConductivity()
+ RANSParentType::eddyThermalConductivity();
}
/*!
* \brief Calculate and set the dynamic eddy viscosity.
* \brief Returns the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$.
*/
void calculateEddyViscosity()
Scalar calculateEddyViscosity()
{
Scalar kinematicEddyViscosity
= cMu() * fMu() * turbulentKineticEnergy() * turbulentKineticEnergy()
/ dissipationTilde();
// if (std::isnan(kinematicEddyViscosity) || dissipationTilde() < 1e-8)
// {
Dune::dinfo /*<< " scv.dofPosition() " << dofPosition_
*/<< " cMu() " << cMu()
<< " fMu() " << fMu()
<< " turbulentKineticEnergy() " << turbulentKineticEnergy()
<< " dissipationTilde() " << dissipationTilde()
<< " kinematicEddyViscosity " << kinematicEddyViscosity
<< std::endl;
// }
ParentType::setKinematicEddyViscosity(kinematicEddyViscosity);
return cMu() * fMu() * turbulentKineticEnergy() * turbulentKineticEnergy()
/ dissipationTilde() * NavierStokesParentType::density();
}
/*!
......@@ -181,15 +182,15 @@ public:
const Scalar reT() const
{
return turbulentKineticEnergy() * turbulentKineticEnergy()
/ ParentType::kinematicViscosity() / dissipationTilde();
/ RANSParentType::kinematicViscosity() / dissipationTilde();
}
//! \brief Returns the \$f Re_\textrm{y} \$f value
const Scalar reY() const
{
using std::sqrt;
return sqrt(turbulentKineticEnergy()) * ParentType::wallDistance()
/ ParentType::kinematicViscosity();
return sqrt(turbulentKineticEnergy()) * RANSParentType::wallDistance()
/ RANSParentType::kinematicViscosity();
}
//! \brief Returns the \$f C_\mu \$f constant
......@@ -255,8 +256,8 @@ public:
else if (lowReKEpsilonModel_ == LowReKEpsilonModels::lamBremhorst)
return 0.0;
else // use default model by Chien
return 2.0 * ParentType::kinematicViscosity() * turbulentKineticEnergy()
/ ParentType::wallDistance() / ParentType::wallDistance();
return 2.0 * RANSParentType::kinematicViscosity() * turbulentKineticEnergy()
/ RANSParentType::wallDistance() / RANSParentType::wallDistance();
}
//! \brief Returns the \$f f_\mu \$f value
......@@ -270,7 +271,7 @@ public:
return pow((1.0 - exp(-0.0165 * reY())), 2.0)
* (1.0 + 20.5 / reT());
else // use default model by Chien
return 1.0 - exp(-0.0115 * ParentType::yPlus());
return 1.0 - exp(-0.0115 * RANSParentType::yPlus());
}
//! \brief Returns the \$f f_1 \$f value
......@@ -306,13 +307,14 @@ public:
else if (lowReKEpsilonModel_ == LowReKEpsilonModels::lamBremhorst)
return 0.0;
else // use default model by Chien
return -2.0 * ParentType::kinematicViscosity() * dissipationTilde()
/ ParentType::wallDistance() / ParentType::wallDistance()
* exp(-0.5 * ParentType::yPlus());
return -2.0 * RANSParentType::kinematicViscosity() * dissipationTilde()
/ RANSParentType::wallDistance() / RANSParentType::wallDistance()
* exp(-0.5 * RANSParentType::yPlus());
}
protected:
int lowReKEpsilonModel_;
Scalar dynamicEddyViscosity_;
Scalar turbulentKineticEnergy_;
Scalar dissipationTilde_;
Scalar storedTurbulentKineticEnergy_;
......@@ -320,64 +322,6 @@ protected:
Scalar stressTensorScalarProduct_;
};
/*!
* \ingroup LowReKEpsilonModel
* \brief Volume variables for the non-isothermal single-phase 0-Eq. model.
*/
template <class TypeTag>
class LowReKEpsilonVolumeVariablesImplementation<TypeTag, true>
: virtual public LowReKEpsilonVolumeVariablesImplementation<TypeTag, false>,
virtual public RANSVolumeVariablesImplementation<TypeTag, true>
{
using ParentTypeNonIsothermal = RANSVolumeVariablesImplementation<TypeTag, true>;
using ParentTypeIsothermal = LowReKEpsilonVolumeVariablesImplementation<TypeTag, false>;
using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
public:
/*!
* \brief Update all quantities for a given control volume
*
* \param elemSol A vector containing all primary variables connected to the element
* \param problem The object specifying the problem which ought to be simulated
* \param element An element which contains part of the control volume
* \param scv The sub-control volume
*/
template<class ElementSolution>
void update(const ElementSolution &elemSol,
const Problem &problem,
const Element &element,
const SubControlVolume &scv)
{
ParentTypeNonIsothermal::update(elemSol, problem, element, scv);
updateRANSProperties(elemSol, problem, element, scv);
}
/*!
* \brief Update all turbulent quantities for a given control volume
*
* Wall and roughness related quantities are stored. Eddy viscosity is set.
*
* \param elemSol A vector containing all primary variables connected to the element
* \param problem The object specifying the problem which ought to be simulated
* \param element An element which contains part of the control volume
* \param scv The sub-control volume
*/
template<class ElementSolution>
void updateRANSProperties(const ElementSolution &elemSol,
const Problem &problem,
const Element &element,
const SubControlVolume& scv)
{
ParentTypeIsothermal::updateRANSProperties(elemSol, problem, element, scv);
ParentTypeNonIsothermal::calculateEddyThermalConductivity();
}
};
}
#endif
......@@ -95,9 +95,6 @@ public:
uStar_ = max(uStar_, 1e-10); // zero values lead to numerical problems in some turbulence models
yPlus_ = wallDistance_ * uStar_ / problem.kinematicViscosity_[elementID_];
uPlus_ = velocity_[flowNormalAxis] / uStar_;
// get the dynamic eddy viscosity from the specific RANS implementation
dynamicEddyViscosity_ = asImp_().calculateEddyViscosity(elemSol, problem, element, scv);
}
/*!
......@@ -153,7 +150,7 @@ public:
* control volume.
*/
Scalar kinematicEddyViscosity() const
{ return dynamicEddyViscosity() / asImp_().density(); }
{ return asImp_().dynamicEddyViscosity() / asImp_().density(); }
/*!
* \brief Return the kinematic viscosity \f$\mathrm{[m^2/s]}\f$ of the fluid within the
......@@ -162,18 +159,10 @@ public:
Scalar kinematicViscosity() const
{ return asImp_().viscosity() / asImp_().density(); }
/*!
* \brief Return the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$ of the flow within the
* control volume.
*/
Scalar dynamicEddyViscosity() const
{ return dynamicEddyViscosity_; }
protected:
DimVector velocity_;
DimVector velocityMaximum_;
DimMatrix velocityGradients_;
Scalar dynamicEddyViscosity_;
std::size_t elementID_;
std::size_t wallElementID_;
Scalar wallDistance_;
......
......@@ -72,18 +72,41 @@ public:
const SubControlVolume& scv)
{
NavierStokesParentType::update(elemSol, problem, element, scv);
updateRANSProperties(elemSol, problem, element, scv);
}
/*!
* \brief Update all turbulent quantities for a given control volume
*
* \param elemSol A vector containing all primary variables connected to the element
* \param problem The object specifying the problem which ought to be simulated
* \param element An element which contains part of the control volume
* \param scv The sub-control volume
*/
template<class ElementSolution, class Problem, class Element, class SubControlVolume>
void updateRANSProperties(const ElementSolution &elemSol,
const Problem &problem,
const Element &element,
const SubControlVolume& scv)
{
RANSParentType::updateRANSProperties(elemSol, problem, element, scv);
additionalRoughnessLength_ = problem.additionalRoughnessLength_[RANSParentType::elementID()];
yPlusRough_ = wallDistanceRough() * RANSParentType::uStar() / RANSParentType::kinematicViscosity();
dynamicEddyViscosity_ = calculateEddyViscosity(elemSol, problem, element, scv);
}
/*!
* \brief Return the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$ of the flow
*/
Scalar dynamicEddyViscosity() const
{ return dynamicEddyViscosity_; }
/*!
* \brief Return the effective dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
* control volume.
*/
Scalar effectiveViscosity() const
{
return NavierStokesParentType::viscosity()
+ RANSParentType::dynamicEddyViscosity();
}
{ return NavierStokesParentType::viscosity() + dynamicEddyViscosity(); }
/*!
* \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
......@@ -110,9 +133,6 @@ public:
const Element &element,
const SubControlVolume& scv)
{
additionalRoughnessLength_ = problem.additionalRoughnessLength_[RANSParentType::elementID()];
yPlusRough_ = wallDistanceRough() * RANSParentType::uStar() / RANSParentType::kinematicViscosity();
using std::abs;
using std::exp;
using std::sqrt;
......@@ -169,7 +189,7 @@ public:
{ return NavierStokesParentType::fluidState_; }
protected:
Scalar dynamicEddyViscosity_;
Scalar additionalRoughnessLength_;
Scalar yPlusRough_;
};
......
......@@ -15,6 +15,6 @@ dune_add_test(NAME test_channel_zeroeq2cni
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_channel_zeroeq2cni.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_channel_zeroeq2cni-00016.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_channel_zeroeq2cni-00019.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_zeroeq2cni test_channel_zeroeq2cni.input")
target_compile_definitions(test_channel_zeroeq2cni PUBLIC "NONISOTHERMAL=1")
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment