freeflow/rans/oneeq/volumevariables.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_ONEEQ_VOLUME_VARIABLES_HH
#define DUMUX_ONEEQ_VOLUME_VARIABLES_HH
 
#include <dune/common/math.hh>
#include <dumux/common/parameters.hh>
#include <dumux/freeflow/rans/volumevariables.hh>
 
namespace Dumux {
 
template <class Traits, class NSVolumeVariables>
class OneEqVolumeVariables
:  public RANSVolumeVariables<Traits, NSVolumeVariables>
{
    using RANSParentType = RANSVolumeVariables<Traits, NSVolumeVariables>;
 
    using Scalar = typename Traits::PrimaryVariables::value_type;
    using DimVector = Dune::FieldVector<Scalar, Traits::ModelTraits::dim()>;
 
public:
    using Indices = typename Traits::ModelTraits::Indices;
 
    template<class ElementSolution, class Problem, class Element, class SubControlVolume>
    void update(const ElementSolution &elemSol,
                const Problem &problem,
                const Element &element,
                const SubControlVolume& scv)
    {
        RANSParentType::updateNavierStokesVolVars(elemSol, problem, element, scv);
        updateRANSProperties(elemSol, problem, element, scv);
    }
 
    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);
        viscosityTilde_ = elemSol[0][Indices::viscosityTildeIdx];
        storedViscosityTilde_ = problem.storedViscosityTilde(RANSParentType::elementIdx());
        storedViscosityTildeGradient_ = problem.storedViscosityTildeGradient(RANSParentType::elementIdx());
        stressTensorScalarProduct_ = problem.stressTensorScalarProduct(RANSParentType::elementIdx());
        vorticityTensorScalarProduct_ = problem.vorticityTensorScalarProduct(RANSParentType::elementIdx());
        if (problem.useStoredEddyViscosity())
            RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity(RANSParentType::elementIdx()));
        else
            RANSParentType::setDynamicEddyViscosity_(calculateEddyViscosity());
        RANSParentType::calculateEddyDiffusivity(problem);
        RANSParentType::calculateEddyThermalConductivity(problem);
    }
 
    Scalar calculateEddyViscosity()
    { return viscosityTilde() * fv1() *  RANSParentType::density(); }
 
    Scalar viscosityTilde() const
    { return viscosityTilde_; }
 
    Scalar storedViscosityTilde() const
    { return storedViscosityTilde_; }
 
    DimVector storedViscosityTildeGradient() const
    { return storedViscosityTildeGradient_; }
 
    Scalar stressTensorScalarProduct() const
    { return stressTensorScalarProduct_; }
 
    Scalar fv1() const
    {
        return viscosityRatio() * viscosityRatio() * viscosityRatio()
               / (viscosityRatio() * viscosityRatio() * viscosityRatio()
                  + cv1() * cv1() * cv1());
    }
 
    Scalar fv2() const
    { return 1.0 - viscosityRatio() / (1.0 + viscosityRatio() * fv1()); }
 
    Scalar ft2() const
    {
        using std::exp;
        // the trip correction term is dropped according to Versteeg2009 and Wilcox2006
        // return ct3() * exp(-ct4() * viscosityRatio()  * viscosityRatio());
        return 0.0;
    }
 
    Scalar fW() const
    {
        using std::pow;
        using Dune::power;
        return g() * pow(1.0 + power(cw3(), 6) / (power(g(), 6) + power(cw3(), 6)), 1.0/6.0);
    }
 
    Scalar g() const
    {
        using Dune::power;
        return r() + cw2() * (power(r(), 6) - r());
    }
 
    Scalar r() const
    {
        using std::min;
        return min(10.0,
                   viscosityTilde() / stressTensorScalarProductTilde()
                   / RANSParentType::karmanConstant() / RANSParentType::karmanConstant()
                   / RANSParentType::wallDistance() / RANSParentType::wallDistance());
    }
 
    Scalar viscosityRatio() const
    { return viscosityTilde() / RANSParentType::kinematicViscosity(); }
 
    /*
     * ! \brief Returns a modified version of the stress tensor scalar product
     *
     * According to <a href="https://turbmodels.larc.nasa.gov/spalart.html">NASA</a>
     * this term should never be zero and different limiters might be used.
     * The Implementation uses the one proposed in:
     * Allmaras, S. R., Johnson, F. T., and Spalart, P. R.,
     * "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model," ICCFD7-1902
     */
    Scalar stressTensorScalarProductTilde() const
    {
        // original form
        // return vorticityMagnitude()
        //        + viscosityTilde() * fv2()
        //          / RANSParentType::karmanConstant() / RANSParentType::karmanConstant()
        //          / RANSParentType::wallDistance() / RANSParentType::wallDistance();
 
        // limiter form, literature source see above
        Scalar sBar = viscosityTilde() * fv2()
                      / RANSParentType::karmanConstant() / RANSParentType::karmanConstant()
                      / RANSParentType::wallDistance() / RANSParentType::wallDistance();
        return sBar < -c2() * vorticityMagnitude()
               ? vorticityMagnitude()
                 + (vorticityMagnitude() * (c2() * c2() * vorticityMagnitude() + c3() * sBar))
                   / ((c3() - 2.0 * c2()) * vorticityMagnitude() - sBar)
               : vorticityMagnitude() + sBar;
    }
 
    Scalar vorticityMagnitude() const
    {
        using std::sqrt;
        return sqrt(2.0 * vorticityTensorScalarProduct_);
    }
 
    Scalar c2() const
    { return 0.7; }
 
    Scalar c3() const
    { return 0.9; }
 
    Scalar sigma() const
    { return 2.0/3.0; }
 
    Scalar cb1() const
    { return 0.1355; }
 
    Scalar cb2() const
    { return 0.622; }
 
    Scalar cv1() const
    { return 7.1; }
 
    Scalar ct3() const
    { return 1.2; }
 
    Scalar ct4() const
    { return 0.5; }
 
    Scalar cw1() const
    {
        return cb1() / RANSParentType::karmanConstant() / RANSParentType::karmanConstant()
               + (1.0 + cb2()) / sigma();
    }
 
    Scalar cw2() const
    { return 0.3; }
 
    Scalar cw3() const
    { return 2.0; }
 
protected:
    Scalar viscosityTilde_ = 0.0;
    Scalar storedViscosityTilde_ = 0.0;
    DimVector storedViscosityTildeGradient_ = DimVector(0.0);
    Scalar stressTensorScalarProduct_;
    Scalar vorticityTensorScalarProduct_;
};
 
} // end namespace Dumux
 
#endif