diff --git a/src/atmosphericModels/Make/files b/src/atmosphericModels/Make/files index b68f5c7b7e..9e4f3018a5 100644 --- a/src/atmosphericModels/Make/files +++ b/src/atmosphericModels/Make/files @@ -1,5 +1,6 @@ /* Models */ -atmosphericTurbulentTransportModels.C +turbulenceModels/atmosphericTurbulentTransportModels.C +turbulenceModels/atmosphericTurbulentFluidThermoModels.C porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C diff --git a/src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C b/src/atmosphericModels/turbulenceModels/RAS/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C similarity index 100% rename from src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C rename to src/atmosphericModels/turbulenceModels/RAS/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C diff --git a/src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.H b/src/atmosphericModels/turbulenceModels/RAS/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.H similarity index 100% rename from src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.H rename to src/atmosphericModels/turbulenceModels/RAS/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.H diff --git a/src/atmosphericModels/turbulenceModels/RAS/kL/kL.C b/src/atmosphericModels/turbulenceModels/RAS/kL/kL.C new file mode 100644 index 0000000000..94e393c2ee --- /dev/null +++ b/src/atmosphericModels/turbulenceModels/RAS/kL/kL.C @@ -0,0 +1,494 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 3 of the License, or + (at your option) any later version. + + OpenFOAM 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 OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "kL.H" +#include "fvOptions.H" +#include "bound.H" +#include "gravityMeshObject.H" +#include "wallDist.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace RASModels +{ + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +template +tmp kL::Cmu() const +{ + // (A:Eq. 31) + return (0.556 + 0.108*Rt_)/(1.0 + 0.308*Rt_ + 0.00837*sqr(Rt_)); +} + + +template +tmp kL::CmuPrime() const +{ + // (A:Eq. 32) + return 0.556/(1.0 + 0.277*Rt_); +} + + +template +tmp kL::nutPrime() const +{ + // (A:Eq. 12) + return CmuPrime()*sqrt(k_)*L_; +} + + +template +tmp kL::epsilonCanopy() const +{ + const auto* CdPtr = + this->mesh_.template findObject("plantCd"); + const auto* LADPtr = + this->mesh_.template findObject("leafAreaDensity"); + const volVectorField& U = this->U_; + + if (CdPtr && LADPtr) + { + const auto& Cd = *CdPtr; + const auto& LAD = *LADPtr; + + // (W:Eq. 13) + return Cd*LAD*mag(U)*k_; + } + + return tmp::New + ( + IOobject + ( + IOobject::groupName("epsilonCanopy", this->alphaRhoPhi_.group()), + this->runTime_.timeName(), + this->mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + this->mesh_, + dimensionedScalar(sqr(dimLength)/pow3(dimTime), Zero) + ); +} + + +template +tmp kL::epsilon() const +{ + // (W:Eq. 13) + tmp tepsilonCanopy = epsilonCanopy(); + + // (A:Eq. 19) + tmp tepsilonPlain = pow3(Cmu0_)*pow(k_, 1.5)/L_; + + // (W:Eq. 13) + tmp tepsilon = max(tepsilonPlain, tepsilonCanopy); + volScalarField& epsilon = tepsilon.ref(); + bound(epsilon, this->epsilonMin_); + + return tepsilon; +} + + +template +tmp kL::canopyHeight() const +{ + const auto* canopyHeightPtr = + this->mesh_.template findObject("canopyHeight"); + + if (canopyHeightPtr) + { + const auto& canopyHeight = *canopyHeightPtr; + return canopyHeight; + } + + return tmp::New + ( + IOobject + ( + IOobject::groupName("canopyHeight", this->alphaRhoPhi_.group()), + this->runTime_.timeName(), + this->mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + this->mesh_, + dimensionedScalar(dimLength, Zero) + ); +} + + +template +tmp kL::L() const +{ + // (A:Eq. 22) + const volScalarField Lplain(kappa_*y_); + + // Increase roughness for canopy (forest, vegetation etc) + tmp tLcanopy = kappa_*canopyHeight(); + const volScalarField& Lcanopy = tLcanopy; + + // (W:Eq. 16) + return max(Lcanopy, Lplain); +} + + +template +void kL::stratification(const volScalarField& fVB) +{ + tmp tLg = L(); + const volScalarField& Lg = tLg.cref(); + + tmp tcanopyHeight = canopyHeight(); + const volScalarField& canopyHeight = tcanopyHeight; + + tmp tLcanopy = kappa_*canopyHeight; + const volScalarField& Lcanopy = tLcanopy; + + const scalar Cmu0 = Cmu0_.value(); + const scalar CbStable = CbStable_.value(); + const scalar CbUnstable = CbUnstable_.value(); + + forAll(L_, i) + { + if (y_[i] > canopyHeight[i]) + { + if (fVB[i] > 0) + { + // (A:Eq. 23) + const scalar Lb = CbStable*sqrt(k_[i])/sqrt(fVB[i]); + + // (A:Eq. 26) + L_[i] = sqrt(sqr(Lg[i]*Lb)/(sqr(Lg[i]) + sqr(Lb))); + } + else + { + // For unstable/neutral boundary layer (A:p. 80) + // Smoothing function for turbulent Richardson + // number to ensure gentle transition into + // the regime of strong convection + Rt_[i] = + min + ( + max(Rt_[i], -1.0), + Rt_[i] - sqr(Rt_[i] + 1.0)/(Rt_[i] - 1.0) + ); + + // (A:Eq. 28) + L_[i] = + Lg[i] + *sqrt(1.0 - pow(Cmu0, 6.0)*pow(CbUnstable, -2.0)*Rt_[i]); + } + } + else + { + L_[i] = Lcanopy[i]; + } + } + + // Limit characteristic length scale + L_ = min(L_, Lmax_); +} + + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +template +void kL::correctNut() +{ + this->nut_ = Cmu()*sqrt(k_)*L_; + this->nut_.correctBoundaryConditions(); + fv::options::New(this->mesh_).correct(this->nut_); + + BasicTurbulenceModel::correctNut(); +} + + +template +tmp kL::kSource() const +{ + return tmp::New + ( + k_, + dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime + ); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +kL::kL +( + const alphaField& alpha, + const rhoField& rho, + const volVectorField& U, + const surfaceScalarField& alphaRhoPhi, + const surfaceScalarField& phi, + const transportModel& transport, + const word& propertiesName, + const word& type +) +: + eddyViscosity> + ( + type, + alpha, + rho, + U, + alphaRhoPhi, + phi, + transport, + propertiesName + ), + + kappa_ + ( + dimensioned::getOrAddToDict + ( + "kappa", + this->coeffDict_, + 0.41 + ) + ), + sigmak_ + ( + dimensioned::getOrAddToDict + ( + "sigmak", + this->coeffDict_, + 1.0 + ) + ), + beta_ + ( + dimensioned::getOrAddToDict + ( + "beta", + this->coeffDict_, + dimless/dimTemperature, + 3.3e-03 + ) + ), + Cmu0_ + ( + dimensioned::getOrAddToDict + ( + "Cmu0", + this->coeffDict_, + 0.556 + ) + ), + Lmax_ + ( + dimensioned::getOrAddToDict + ( + "Lmax", + this->coeffDict_, + dimLength, + GREAT + ) + ), + CbStable_ + ( + dimensioned::getOrAddToDict + ( + "CbStable", + this->coeffDict_, + 0.25 + ) + ), + CbUnstable_ + ( + dimensioned::getOrAddToDict + ( + "CbUnstable", + this->coeffDict_, + 0.35 + ) + ), + + k_ + ( + IOobject + ( + IOobject::groupName("k", alphaRhoPhi.group()), + this->runTime_.timeName(), + this->mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + this->mesh_ + ), + L_ + ( + IOobject + ( + IOobject::groupName("L", alphaRhoPhi.group()), + this->runTime_.timeName(), + this->mesh_, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + this->mesh_, + dimensionedScalar(dimLength, scalar(1)) + ), + Rt_ + ( + IOobject + ( + IOobject::groupName("Rt", alphaRhoPhi.group()), + this->runTime_.timeName(), + this->mesh_, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + this->mesh_, + dimensionedScalar(dimless, Zero) + ), + g_(meshObjects::gravity::New(this->mesh_.time())), + y_(wallDist::New(this->mesh_).y()) +{ + bound(k_, this->kMin_); + + if (type == typeName) + { + this->printCoeffs(type); + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +bool kL::read() +{ + if (eddyViscosity>::read()) + { + kappa_.readIfPresent(this->coeffDict()); + sigmak_.readIfPresent(this->coeffDict()); + beta_.readIfPresent(this->coeffDict()); + Cmu0_.readIfPresent(this->coeffDict()); + Lmax_.readIfPresent(this->coeffDict()); + CbStable_.readIfPresent(this->coeffDict()); + CbUnstable_.readIfPresent(this->coeffDict()); + + return true; + } + + return false; +} + + +template +void kL::correct() +{ + if (!this->turbulence_) + { + return; + } + + // Construct local convenience references + const alphaField& alpha = this->alpha_; + const rhoField& rho = this->rho_; + const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_; + const volVectorField& U = this->U_; + const volScalarField& nut = this->nut_; + + fv::options& fvOptions(fv::options::New(this->mesh_)); + + eddyViscosity>::correct(); + + // Turbulent kinetic energy production rate + tmp tgradU = fvc::grad(U); + const volScalarField::Internal G + ( + this->GName(), + nut.v()*2*magSqr(dev(symm(tgradU.cref().v()))) + ); + tgradU.clear(); + + // Square of Brunt-Vaisala (buoyancy) frequency + const auto& T = U.mesh().lookupObject("T"); + tmp tfBV = -beta_*(fvc::grad(T) & g_); + const volScalarField& fBV = tfBV.cref(); + + // Sink or source of TKE depending on stratification type (A:Eq. 15) + tmp tPb = -fBV*nutPrime(); + const volScalarField& Pb = tPb.cref(); + + // Turbulent kinetic energy dissipation rate due to plains and canopy + tmp tepsilon = epsilon(); + const volScalarField& epsilon = tepsilon.cref(); + + // Divergence of velocity + tmp tdivU = fvc::div(fvc::absolute(this->phi(), U)); + const volScalarField::Internal& divU = tdivU.cref().v(); + + // Turbulent kinetic energy equation + tmp kEqn + ( + fvm::ddt(alpha, rho, k_) + + fvm::div(alphaRhoPhi, k_) + - fvm::laplacian(alpha*rho*DkEff(), k_) + == + alpha()*rho()*G + + fvm::SuSp((Pb - epsilon)/k_, k_) + - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_) + + kSource() + + fvOptions(alpha, rho, k_) + ); + + tdivU.clear(); + tPb.clear(); + + kEqn.ref().relax(); + fvOptions.constrain(kEqn.ref()); + solve(kEqn); + fvOptions.correct(k_); + bound(k_, this->kMin_); + + // Turbulent Richardson number (A:Eq. 29) + Rt_ = fBV*sqr(k_/tepsilon); + + stratification(fBV); + tfBV.clear(); + + correctNut(); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/atmosphericModels/turbulenceModels/RAS/kL/kL.H b/src/atmosphericModels/turbulenceModels/RAS/kL/kL.H new file mode 100644 index 0000000000..870b38ac69 --- /dev/null +++ b/src/atmosphericModels/turbulenceModels/RAS/kL/kL.H @@ -0,0 +1,313 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 3 of the License, or + (at your option) any later version. + + OpenFOAM 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 OpenFOAM. If not, see . + +Class + Foam::RASModels::kL + +Group + grpRASTurbulence + +Description + A one-equation (turbulent kinetic energy \c k) turbulence closure + model for incompressible and compressible geophysical applications. + + Turbulent kinetic energy (\c k) is computed with a transport equation + and the turbulent length scale (\c L) is computed with an algebraic + expression which depends on the local stratification. + + References: + \verbatim + Standard model (tag:A): + Axell, L. B., & Liungman, O. (2001). + A one-equation turbulence model for geophysical applications: + comparison with data and the k−ε model. + Environmental Fluid Mechanics, 1(1), 71-106. + DOI:10.1023/A:1011560202388 + + Canopy-related models (tag:W): + Wilson, J. D., & Flesch, T. K. (1999). + Wind and remnant tree sway in forest cutblocks. + III. A windflow model to diagnose spatial variation. + Agricultural and Forest Meteorology, 93(4), 259-282. + DOI:10.1016/S0168-1923(98)00121-X + \endverbatim + +Usage + Example by using \c constant/turbulenceProperties: + \verbatim + RAS + { + // Mandatory entries + RASModel kL; + + // Optional entries + kLCoeffs + { + kappa ; + sigmak ; + beta ; + Cmu0 ; + Lmax ; + CbStable ; + CbUnstable ; + } + + // Inherited entries + ... + } + \endverbatim + + where the entries mean: + \table + Property | Description | Type | Reqd | Deflt + RASModel | Type name: kL | word | yes | - + kappa | von Karman constant | scalar | no | 0.41 + sigmak | Empirical model coefficient | scalar | no | 1.0 + beta | Thermal expansion coefficient [1/K] | scalar | no | 3.3e-3 + Cmu0 | Empirical model coefficient | scalar | no | 0.556 + Lmax | Maximum mixing-length scale [m] | scalar | no | GREAT + CbStable | Stable stratification constant | scalar | no | 0.25 + CbUnstable | Unstable stratification constant | scalar | no | 0.35 + \endtable + + The inherited entries are elaborated in: + - \link eddyViscosity.H \endlink + + \heading Input fields (mandatory) + \plaintable + k | Turbulent kinetic energy [m2/s2] + T | Potential temperature [K] + \endplaintable + + \heading Input fields (optional) + \plaintable + canopyHeight | Canopy height [m] + plantCd | Plant canopy drag coefficient [-] + leafAreaDensity | Leaf area density [1/m] + Rt | Turbulent Richardson number [-] + L | Characteristic length scale [m] + \endplaintable + +Note + - Optional input fields can/should be input + by using \c readFields function object. + +SourceFiles + kL.C + +\*---------------------------------------------------------------------------*/ + +#ifndef kL_H +#define kL_H + +#include "RASModel.H" +#include "eddyViscosity.H" +#include "uniformDimensionedFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace RASModels +{ + +/*---------------------------------------------------------------------------*\ + Class kL Declaration +\*---------------------------------------------------------------------------*/ + +template +class kL +: + public eddyViscosity> +{ + // Private Member Functions + + //- Return modified stability function of Launder, + //- calibrated with the results of Hogstrom (A:Eq. 31) + tmp Cmu() const; + + //- Return modified stability function of Launder, + //- calibrated with the results of Hogstrom (A:Eq. 32) + tmp CmuPrime() const; + + //- Return eddy diffusivity (A:Eq. 12) + tmp nutPrime() const; + + //- Return turbulent kinetic energy dissipation rate due to canopy + tmp epsilonCanopy() const; + + //- Return turbulent kinetic energy + //- dissipation rate due to plains and canopy + tmp epsilon() const; + + //- Return canopy height + tmp canopyHeight() const; + + //- Return characteristic length scale + tmp L() const; + + //- Modify characteristic length scale + //- according to the specified stratification + void stratification(const volScalarField& fVB); + + // Generated Methods + + //- No copy construct + kL(const kL&) = delete; + + //- No copy assignment + void operator=(const kL&) = delete; + + +protected: + + // Protected Data + + // Model coefficients + + //- von Karman constant + dimensionedScalar kappa_; + + //- Empirical model coefficient + dimensionedScalar sigmak_; + + //- Thermal expansion coefficient [1/K] + dimensionedScalar beta_; + + //- Empirical model coefficient + dimensionedScalar Cmu0_; + + //- Maximum mixing-length scalar [m] + dimensionedScalar Lmax_; + + //- Stable stratification constant + dimensionedScalar CbStable_; + + //- Unstable stratification constant + dimensionedScalar CbUnstable_; + + + // Fields + + //- Turbulent kinetic energy [m2/s2] + volScalarField k_; + + //- Characteristic length scale [m] + volScalarField L_; + + //- Turbulent Richardson number [-] + volScalarField Rt_; + + //- Gravitational acceleration [m2/s2] + const uniformDimensionedVectorField& g_; + + //- Wall distance + // Note: different to wall distance in parent RASModel + // which is for near-wall cells only + const volScalarField& y_; + + + // Protected Member Functions + + //- Correct the turbulence viscosity + virtual void correctNut(); + + //- Add explicit source for turbulent kinetic energy + virtual tmp kSource() const; + + +public: + + typedef typename BasicTurbulenceModel::alphaField alphaField; + typedef typename BasicTurbulenceModel::rhoField rhoField; + typedef typename BasicTurbulenceModel::transportModel transportModel; + + + //- Runtime type information + TypeName("kL"); + + + // Constructors + + //- Construct from components + kL + ( + const alphaField& alpha, + const rhoField& rho, + const volVectorField& U, + const surfaceScalarField& alphaRhoPhi, + const surfaceScalarField& phi, + const transportModel& transport, + const word& propertiesName = turbulenceModel::propertiesName, + const word& type = typeName + ); + + + //- Destructor + virtual ~kL() = default; + + + // Member Functions + + //- Re-read model coefficients if they have changed + virtual bool read(); + + //- Return the effective diffusivity for k + tmp DkEff() const + { + return tmp::New + ( + "DkEff", + (this->nut_/sigmak_ + this->nu()) + ); + } + + //- Return the turbulent kinetic energy field + virtual tmp k() const + { + return k_; + } + + //- Solve the turbulence equations and correct the turbulence viscosity + virtual void correct(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "kL.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/atmosphericModels/turbulenceModels/atmosphericTurbulentFluidThermoModels.C b/src/atmosphericModels/turbulenceModels/atmosphericTurbulentFluidThermoModels.C new file mode 100644 index 0000000000..d813a5dbb1 --- /dev/null +++ b/src/atmosphericModels/turbulenceModels/atmosphericTurbulentFluidThermoModels.C @@ -0,0 +1,40 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 3 of the License, or + (at your option) any later version. + + OpenFOAM 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 OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "turbulentFluidThermoModels.H" + +// -------------------------------------------------------------------------- // +// RAS models +// -------------------------------------------------------------------------- // + +#include "kEpsilonLopesdaCosta.H" +makeRASModel(kEpsilonLopesdaCosta); + +#include "kL.H" +makeRASModel(kL); + +// ************************************************************************* // diff --git a/src/atmosphericModels/atmosphericTurbulentTransportModels.C b/src/atmosphericModels/turbulenceModels/atmosphericTurbulentTransportModels.C similarity index 95% rename from src/atmosphericModels/atmosphericTurbulentTransportModels.C rename to src/atmosphericModels/turbulenceModels/atmosphericTurbulentTransportModels.C index 2f1373b65c..48e32541c9 100644 --- a/src/atmosphericModels/atmosphericTurbulentTransportModels.C +++ b/src/atmosphericModels/turbulenceModels/atmosphericTurbulentTransportModels.C @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2018 OpenFOAM Foundation + Copyright (C) 2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -34,4 +35,7 @@ License #include "kEpsilonLopesdaCosta.H" makeRASModel(kEpsilonLopesdaCosta); +#include "kL.H" +makeRASModel(kL); + // ************************************************************************* //