freeflow/rans/twoeq/sst/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_SST_VOLUME_VARIABLES_HH
#define DUMUX_SST_VOLUME_VARIABLES_HH
 
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/freeflow/turbulencemodel.hh>
#include <dumux/freeflow/rans/volumevariables.hh>
 
namespace Dumux {
 
template <class Traits, class NSVolumeVariables>
class SSTVolumeVariables
:  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);
 
        turbulentKineticEnergy_ = elemSol[0][Indices::turbulentKineticEnergyIdx];
        dissipation_ = elemSol[0][Indices::dissipationIdx];
 
        storedDissipation_ = problem.storedDissipation(RANSParentType::elementIdx());
        storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy(RANSParentType::elementIdx());
 
        storedDissipationGradient_ = problem.storedDissipationGradient(RANSParentType::elementIdx());
        storedTurbulentKineticEnergyGradient_ = problem.storedTurbulentKineticEnergyGradient(RANSParentType::elementIdx());
 
        stressTensorScalarProduct_ = problem.stressTensorScalarProduct(RANSParentType::elementIdx());
        vorticityTensorScalarProduct_ = problem.vorticityTensorScalarProduct(RANSParentType::elementIdx());
        wallDistance_ = problem.wallDistance(RANSParentType::elementIdx());
        kinematicViscosity_ = problem.kinematicViscosity(RANSParentType::elementIdx());
 
        if (problem.useStoredEddyViscosity())
            RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity(RANSParentType::elementIdx()));
        else
            RANSParentType::setDynamicEddyViscosity_(calculateEddyViscosity(problem));
 
        RANSParentType::calculateEddyDiffusivity(problem);
        RANSParentType::calculateEddyThermalConductivity(problem);
    }
 
    template<class Problem>
    Scalar calculateEddyViscosity(const Problem& problem)
    {
        using std::sqrt;
        using std::max;
 
        if(problem.sstModelVersion() == SSTModel::BSL)
            return turbulentKineticEnergy() / dissipation() * RANSParentType::density();
        else if(problem.sstModelVersion() == SSTModel::SST)
        {
            const Scalar dividend = a1SST() * turbulentKineticEnergy();
            const Scalar possibleMax1 = a1SST()*dissipation();
            const Scalar possibleMax2 = absoluteValueVorticity() * F2();
            const Scalar divisor = std::max(possibleMax1,possibleMax2);
            return dividend / divisor * RANSParentType::density();
        }
        else
            DUNE_THROW(Dune::NotImplemented, "\nThis SST Model is not implemented.\n");
    }
 
 
    Scalar turbulentKineticEnergy() const
    { return turbulentKineticEnergy_; }
 
    Scalar dissipation() const
    { return dissipation_; }
 
    Scalar storedTurbulentKineticEnergy() const
    { return storedTurbulentKineticEnergy_; }
 
    Scalar storedDissipation() const
    { return storedDissipation_; }
 
    DimVector storedTurbulentKineticEnergyGradient() const
    { return storedTurbulentKineticEnergyGradient_; }
 
    DimVector storedDissipationGradient() const
    { return storedDissipationGradient_; }
 
    Scalar stressTensorScalarProduct() const
    { return stressTensorScalarProduct_; }
 
    Scalar vorticityTensorScalarProduct() const
    { return vorticityTensorScalarProduct_; }
 
    Scalar kinematicViscosity() const
    { return kinematicViscosity_; }
 
    Scalar wallDistance() const
    { return wallDistance_; }
 
    const Scalar betaOmega() const
    { return 0.0708; }
 
    const Scalar alpha() const
    { return 0.52; }
 
    const Scalar sigmaK() const
    { return 0.6; }
 
    const Scalar sigmaOmega() const
    { return 0.5; }
 
    const Scalar betaK() const
    { return 0.09; }
 
    const Scalar sigmaK1BSL() const
    { return 0.5; }
 
    const Scalar sigmaK1SST() const
    { return 0.85; }
 
    const Scalar sigmaK2() const
    { return 1.0; }
 
    const Scalar sigmaOmega1BSL() const
    { return 0.5; }
 
    const Scalar sigmaOmega1SST() const
    { return 0.5; }
 
    const Scalar sigmaOmega2() const
    { return 0.856; }
 
    const Scalar beta1BSL() const
    { return 0.0750; }
 
    const Scalar beta1SST() const
    { return 0.0750; }
 
    const Scalar beta2() const
    { return 0.0820; }
 
    const Scalar betaStar1BSL() const
    { return 0.09; }
 
    const Scalar betaStar1SST() const
    { return 0.09; }
 
    const Scalar betaStar2() const
    { return 0.09; }
 
    const Scalar kappa1BSL() const
    { return 0.41; }
 
    const Scalar kappa1SST() const
    { return 0.41; }
 
    const Scalar kappa2() const
    { return 0.41; }
 
    const Scalar gamma1BSL() const
    { return (beta1BSL()/betaStar1BSL()) - ((sigmaOmega1BSL()*kappa1BSL()*kappa1BSL())/(std::sqrt(betaStar1BSL()))); }
 
    const Scalar gamma1SST() const
    { return (beta1SST()/betaStar1SST()) - ((sigmaOmega1SST()*kappa1SST()*kappa1SST())/(std::sqrt(betaStar1SST()))); }
 
    const Scalar gamma2() const
    { return (beta2()/betaStar2()) - ((sigmaOmega2()*kappa2()*kappa2())/(std::sqrt(betaStar2()))); }
 
    const Scalar a1SST() const
    { return 0.31; }
 
    Scalar absoluteValueVorticity() const
    { return std::sqrt(2.0 * vorticityTensorScalarProduct()); }
 
    Scalar F1() const
    {
        Scalar gradientProduct = 0.0;
        for (unsigned int i = 0; i < Traits::ModelTraits::dim(); ++i){
            gradientProduct += storedTurbulentKineticEnergyGradient()[i] * storedDissipationGradient()[i];}
 
        Scalar positiveCrossDiffusion = 2.0 * RANSParentType::density() * sigmaOmega2() / storedDissipation() * gradientProduct;
        positiveCrossDiffusion = std::max(positiveCrossDiffusion, 1e-20);
 
        const Scalar possibleMin2 = (4.0 * RANSParentType::density() * sigmaOmega2() * storedTurbulentKineticEnergy())
                                  / (positiveCrossDiffusion * wallDistance() * wallDistance());
 
        const Scalar possibleMax1 = (std::sqrt(storedTurbulentKineticEnergy())) / (0.09 * storedDissipation() * wallDistance());
        const Scalar possibleMax2 = (500.00 * kinematicViscosity()) / (wallDistance() * wallDistance() * storedDissipation());
        const Scalar possibleMin1 = std::max(possibleMax1, possibleMax2);
 
        const Scalar argument = std::min(possibleMin1,possibleMin2);
 
        return tanh(argument * argument * argument * argument);
    }
 
 
    Scalar F2() const
    {
        const Scalar possibleMax1 = 2.0 * (std::sqrt(storedTurbulentKineticEnergy())) / (0.09 * storedDissipation() * wallDistance());
        const Scalar possibleMax2 = (500.0 * kinematicViscosity()) / (wallDistance() * wallDistance() * storedDissipation());
        const Scalar argument = std::max(possibleMax1, possibleMax2);
 
        return tanh(argument * argument);
    }
 
    //------------------------Constants for SST--------------------------------
    const Scalar sigmaKSST() const
    { return F1()*sigmaK1SST() + (1-F1())*sigmaK2(); }
 
    const Scalar sigmaOmegaSST() const
    { return F1()*sigmaOmega1SST() + (1-F1())*sigmaOmega2(); }
 
    const Scalar betaSST() const
    { return F1()*beta1SST() + (1-F1())*beta2(); }
 
    const Scalar betaStarSST() const
    { return F1()*betaStar1SST() + (1-F1())*betaStar2(); }
 
    const Scalar kappaSST() const
    { return F1()*kappa1SST() + (1-F1())*kappa2(); }
 
    const Scalar gammaSST() const
    { return F1()*gamma1SST() + (1-F1())*gamma2(); }
 
 
 
    //------------------------Constants for BSL--------------------------------
    const Scalar sigmaKBSL() const
    { return F1()*sigmaK1BSL() + (1-F1())*sigmaK2(); }
 
    const Scalar sigmaOmegaBSL() const
    { return F1()*sigmaOmega1BSL() + (1-F1())*sigmaOmega2(); }
 
    const Scalar betaBSL() const
    { return F1()*beta1BSL() + (1-F1())*beta2(); }
 
    const Scalar betaStarBSL() const
    { return F1()*betaStar1BSL() + (1-F1())*betaStar2(); }
 
    const Scalar kappaBSL() const
    { return F1()*kappa1BSL() + (1-F1())*kappa2(); }
 
    const Scalar gammaBSL() const
    { return F1()*gamma1BSL() + (1-F1())*gamma2(); }
 
 
protected:
    Scalar betaOmega_ = 0.0;
 
    Scalar dissipation_ = 0.0;
    Scalar storedDissipation_ = 0.0;
    DimVector storedDissipationGradient_ = DimVector(0.0);
    Scalar turbulentKineticEnergy_ = 0.0;
    Scalar storedTurbulentKineticEnergy_ = 0.0;
    DimVector storedTurbulentKineticEnergyGradient_ = DimVector(0.0);
 
    Scalar stressTensorScalarProduct_ = 0.0;
    Scalar vorticityTensorScalarProduct_ = 0.0;
    Scalar wallDistance_ = 0.0;
    Scalar kinematicViscosity_ = 0.0;
};
}
 
#endif