Skip to content
Snippets Groups Projects
Commit 2a4fedcf authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

Merge branch 'feature/rans-include-density-changes' into 'master' 

[rans] add density in all terms of the RANS models

See merge request !2139
parents 85ad1b9a 70d50b23
No related branches found
No related tags found
1 merge request!2139[rans] add density in all terms of the RANS models
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 823 additions and 803 deletions
CHANGELOG.md
+ 1
− 0
View file @ 2a4fedcf
......@@ -18,6 +18,7 @@ Differences Between DuMu<sup>x</sup> 3.3 and DuMu<sup>x</sup> 3.2
- The usage of laws for pc-Sw and kr-Sw has been completely revised. A caller does not have to pass a `parameters` object to the laws anymore. The `spatialParams` now provide a `fluidMatrixInteraction` function which bundles an arbitrary number of
different interaction laws such as a pc-Sw and kr-Sw curve and interfacial areas.
New pre-cached spline laws were added which can help to increase efficiency. The usage of the old interface is deprecated and warnings will be raised. The old interface will be removed after the release of 3.3.
- The RANS models now include variable densities. Compositional or nonisothermal RANS models could produce slightly different, more accurate, results.
### Immediate interface changes not allowing/requiring a deprecation period:
- __Flash/Constraintsolver__: The flashes depending on material laws are immediately required to use new-style material laws (fluidMatrixInteraction interface in spatialparams)
......
dumux/freeflow/rans/oneeq/model.hh
+ 6
− 6
View file @ 2a4fedcf
......@@ -36,12 +36,12 @@
* term which account for the transition or trip, is dropped from the original equations,
* such that the balance equation simplifies to:
* \f[
* \frac{\partial \tilde{\nu}}{\partial t}
* + \nabla \cdot \left( \tilde{\nu} \textbf{v} \right)
* - c_\text{b1} \tilde{S} \tilde{\nu}
* - \frac{1}{\sigma_{\tilde{\nu}}} \nabla \cdot \left( \left[ \nu + \tilde{\nu} \right] \nabla \tilde{\nu} \right)
* - \frac{c_\text{b2}}{\sigma_{\tilde{\nu}}} \left| \nabla \tilde{\nu} \right|^2
* + c_\text{w1} f_\text{w} \frac{\tilde{\nu}^2}{y^2}
* \frac{\partial \tilde{\nu}\varrho}{\partial t}
* + \nabla \cdot \left( \tilde{\nu} \varrho \textbf{v} \right)
* - c_\text{b1} \tilde{S} \tilde{\nu} \varrho
* - \frac{1}{\sigma_{\tilde{\nu}}} \nabla \cdot \left( \left[ \mu + \tilde{\nu} \varrho \right] \nabla \tilde{\nu} \right)
* - \frac{c_\text{b2}}{\sigma_{\tilde{\nu}}} \varrho \left| \nabla \tilde{\nu} \right|^2
* + c_\text{w1} f_\text{w} \varrho \frac{\tilde{\nu}^2}{y^2}
* = 0
* \f]
*
......
dumux/freeflow/rans/oneeq/staggered/fluxvariables.hh
+ 3
− 3
View file @ 2a4fedcf
......@@ -96,7 +96,7 @@ public:
// calculate advective flux
auto upwindTermK = [](const auto& volVars)
{
return volVars.viscosityTilde();
return volVars.viscosityTilde() * volVars.density();
};
flux[viscosityTildeEqIdx]
......@@ -109,9 +109,9 @@ public:
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
// effective diffusion coefficients
Scalar insideCoeff = (insideVolVars.kinematicViscosity() + insideVolVars.viscosityTilde())
Scalar insideCoeff = (insideVolVars.viscosity() + insideVolVars.viscosityTilde() * insideVolVars.density())
/ insideVolVars.sigma();
Scalar outsideCoeff = (outsideVolVars.kinematicViscosity() + outsideVolVars.viscosityTilde())
Scalar outsideCoeff = (outsideVolVars.viscosity() + outsideVolVars.viscosityTilde() * outsideVolVars.density())
/ outsideVolVars.sigma();
// scale by extrusion factor
......
dumux/freeflow/rans/oneeq/staggered/localresidual.hh
+ 5
− 4
View file @ 2a4fedcf
......@@ -81,7 +81,7 @@ public:
const VolumeVariables& volVars) const
{
CellCenterPrimaryVariables storage = ParentType::computeStorageForCellCenter(problem, scv, volVars);
storage[viscosityTildeEqIdx] = volVars.viscosityTilde();
storage[viscosityTildeEqIdx] = volVars.viscosityTilde() * volVars.density();
return storage;
}
......@@ -99,18 +99,19 @@ public:
source[viscosityTildeEqIdx] += volVars.cb1() * (1.0 - volVars.ft2())
* volVars.stressTensorScalarProductTilde()
* volVars.viscosityTilde();
* volVars.viscosityTilde() * volVars.density();
source[viscosityTildeEqIdx] -= (volVars.cw1() * volVars.fW()
- volVars.cb1() * volVars.ft2() / problem.karmanConstant() / problem.karmanConstant())
* volVars.viscosityTilde() * volVars.viscosityTilde()
/ volVars.wallDistance() / volVars.wallDistance();
/ volVars.wallDistance() / volVars.wallDistance() * volVars.density();;
for (unsigned int dimIdx = 0; dimIdx < ModelTraits::dim(); ++dimIdx)
{
source[viscosityTildeEqIdx] += volVars.cb2() / volVars.sigma()
* volVars.storedViscosityTildeGradient()[dimIdx]
* volVars.storedViscosityTildeGradient()[dimIdx];
* volVars.storedViscosityTildeGradient()[dimIdx]
* volVars.density();
}
return source;
......
dumux/freeflow/rans/problem.hh
+ 15
− 6
View file @ 2a4fedcf
......@@ -126,7 +126,8 @@ public:
vorticityTensorScalarProduct_.resize(this->gridGeometry().elementMapper().size(), 0.0);
flowDirectionAxis_.resize(this->gridGeometry().elementMapper().size(), fixedFlowDirectionAxis_);
wallNormalAxis_.resize(this->gridGeometry().elementMapper().size(), fixedWallNormalAxis_);
kinematicViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
storedViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
storedDensity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
if ( !(hasParamInGroup(this->paramGroup(), "RANS.IsFlatWallBounded")))
{
......@@ -254,7 +255,7 @@ public:
calculateMaxMinVelocities_();
calculateStressTensor_();
calculateVorticityTensor_();
storeKinematicViscosity_(curSol);
storeViscosities_(curSol);
}
/*!
......@@ -412,8 +413,14 @@ public:
Scalar vorticityTensorScalarProduct(const int elementIdx) const
{ return vorticityTensorScalarProduct_[elementIdx]; }
Scalar storedViscosity(const int elementIdx) const
{ return storedViscosity_[elementIdx]; }
Scalar storedDensity(const int elementIdx) const
{ return storedDensity_[elementIdx]; }
Scalar kinematicViscosity(const int elementIdx) const
{ return kinematicViscosity_[elementIdx]; }
{ return storedViscosity(elementIdx) / storedDensity(elementIdx); }
bool calledUpdateStaticWallProperties = false;
......@@ -675,7 +682,7 @@ private:
}
}
void storeKinematicViscosity_(const SolutionVector& curSol)
void storeViscosities_(const SolutionVector& curSol)
{
// calculate or call all secondary variables
for (const auto& element : elements(this->gridGeometry().gridView()))
......@@ -693,7 +700,8 @@ private:
VolumeVariables volVars;
volVars.update(elemSol, asImp_(), element, scv);
kinematicViscosity_[elementIdx] = volVars.viscosity() / volVars.density();
storedDensity_[elementIdx] = volVars.density();
storedViscosity_[elementIdx] = volVars.viscosity();
}
}
}
......@@ -716,7 +724,8 @@ private:
std::vector<Scalar> stressTensorScalarProduct_;
std::vector<Scalar> vorticityTensorScalarProduct_;
std::vector<Scalar> kinematicViscosity_;
std::vector<Scalar> storedDensity_;
std::vector<Scalar> storedViscosity_;
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
......
dumux/freeflow/rans/twoeq/kepsilon/model.hh
+ 12
− 12
View file @ 2a4fedcf
......@@ -31,27 +31,27 @@
*
* The turbulent kinetic energy balance is:
* \f[
* \frac{\partial \left( k \right)}{\partial t}
* + \nabla \cdot \left( \textbf{v} k \right)
* - \nabla \cdot \left( \left( \nu + \frac{\nu_\text{t}}{\sigma_\text{k}} \right) \nabla k \right)
* - 2 \nu_\text{t} \textbf{S} \cdot \textbf{S}
* + \varepsilon
* \frac{\partial \left( \varrho k \right)}{\partial t}
* + \nabla \cdot \left( \textbf{v} \varhho k \right)
* - \nabla \cdot \left( \left( \mu + \frac{\mu_\text{t}}{\sigma_\text{k}} \right) \nabla k \right)
* - 2 \mu_\text{t} \textbf{S} \cdot \textbf{S}
* + \varrho \varepsilon
* = 0
* \f].
*
* The dissipation balance is:
* \f[
* \frac{\partial \left( \varepsilon \right)}{\partial t}
* + \nabla \cdot \left( \textbf{v} \varepsilon \right)
* - \nabla \cdot \left( \left( \nu + \frac{\nu_\text{t}}{\sigma_{\varepsilon}} \right) \nabla \varepsilon \right)
* - C_{1\varepsilon} \frac{\varepsilon}{k} 2 \nu_\text{t} \textbf{S} \cdot \textbf{S}
* + C_{2\varepsilon} \frac{\varepsilon^2}{k}
* \frac{\partial \left( \varrho \varepsilon \right)}{\partial t}
* + \nabla \cdot \left( \textbf{v} \varrho \varepsilon \right)
* - \nabla \cdot \left( \left( \mu + \frac{\mu_\text{t}}{\sigma_{\varepsilon}} \right) \nabla \varepsilon \right)
* - C_{1\varepsilon} \frac{\varepsilon}{k} 2 \mu_\text{t} \textbf{S} \cdot \textbf{S}
* + C_{2\varepsilon} \varrho \frac{\varepsilon^2}{k}
* = 0
* \f].
*
* The kinematic eddy viscosity \f$ \nu_\text{t} \f$ is:
* The dynamic eddy viscosity \f$ \mu_\text{t} \f$ is:
* \f[
* \nu_\text{t} = C_\mu \frac{k^2}{\tilde{\varepsilon}}
* \mu_\text{t} = \varrho C_\mu \frac{k^2}{\tilde{\varepsilon}}
* \f].
*
* Finally, the model is closed with the following constants:
......
dumux/freeflow/rans/twoeq/kepsilon/problem.hh
+ 3
− 10
View file @ 2a4fedcf
......@@ -102,7 +102,6 @@ public:
ParentType::updateStaticWallProperties();
// update size and initial values of the global vectors
matchingPointIdx_.resize(this->gridGeometry().elementMapper().size(), 0);
storedDensity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
storedDissipation_.resize(this->gridGeometry().elementMapper().size(), 0.0);
storedTurbulentKineticEnergy_.resize(this->gridGeometry().elementMapper().size(), 0.0);
storedDynamicEddyViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
......@@ -138,7 +137,6 @@ public:
VolumeVariables volVars;
volVars.update(elemSol, asImp_(), element, scv);
storedDynamicEddyViscosity_[elementIdx] = volVars.calculateEddyViscosity();
storedDensity_[elementIdx] = volVars.density();
}
}
......@@ -255,7 +253,7 @@ public:
unsigned int wallNormalAxis = asImp_().wallNormalAxis(elementIdx);
unsigned int flowDirectionAxis = asImp_().flowDirectionAxis(elementIdx);
Scalar velocityGradient = asImp_().velocityGradient(elementIdx, flowDirectionAxis, wallNormalAxis);
return mixingLength * mixingLength * abs(velocityGradient) * storedDensity(elementIdx);
return mixingLength * mixingLength * abs(velocityGradient) * asImp_().storedDensity(elementIdx);
}
//! \brief Returns the wall shear stress velocity
......@@ -302,8 +300,7 @@ public:
{
using std::log;
Scalar velocityNominal = uStarNominal(elementIdx) * (1.0 / asImp_().karmanConstant() * log(yPlusNominal(elementIdx)) + 5.0);
return uStarNominal(elementIdx) * uStarNominal(elementIdx)
* velocity / velocityNominal;
return uStarNominal(elementIdx) * uStarNominal(elementIdx) * velocity / velocityNominal;
}
//! \brief Checks whether a wall function should be used
......@@ -328,7 +325,7 @@ public:
{
unsigned int elementIdx = asImp_().gridGeometry().elementMapper().index(element);
return FacePrimaryVariables(asImp_().tangentialMomentumWallFunction(elementIdx, elemFaceVars[scvf].velocitySelf())
* elemVolVars[scvf.insideScvIdx()].density());
* asImp_().storedDensity(elementIdx) );
}
//! \brief Returns the flux for non-isothermal and compositional RANS models
......@@ -469,9 +466,6 @@ public:
return useStoredEddyViscosity;
}
Scalar storedDensity(const int elementIdx) const
{ return storedDensity_[elementIdx]; }
Scalar storedDissipation(const int elementIdx) const
{ return storedDissipation_[elementIdx]; }
......@@ -489,7 +483,6 @@ public:
private:
std::vector<unsigned int> matchingPointIdx_;
std::vector<Scalar> storedDensity_;
std::vector<Scalar> storedDissipation_;
std::vector<Scalar> storedTurbulentKineticEnergy_;
std::vector<Scalar> storedDynamicEddyViscosity_;
......
dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
+ 6
− 15
View file @ 2a4fedcf
......@@ -97,11 +97,11 @@ public:
// calculate advective flux
auto upwindTermK = [](const auto& volVars)
{
return volVars.turbulentKineticEnergy();
return volVars.turbulentKineticEnergy() * volVars.density();
};
auto upwindTermEpsilon = [](const auto& volVars)
{
return volVars.dissipation();
return volVars.dissipation() * volVars.density();
};
flux[turbulentKineticEnergyEqIdx]
......@@ -116,19 +116,10 @@ public:
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
// effective diffusion coefficients
Scalar insideCoeff_k = insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaK();
Scalar outsideCoeff_k = outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaK();
Scalar insideCoeff_e = insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaEpsilon();
Scalar outsideCoeff_e = outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaEpsilon();
static const auto kEpsilonEnableKinematicViscosity_
= getParamFromGroup<bool>(problem.paramGroup(), "KEpsilon.EnableKinematicViscosity", true);
if (kEpsilonEnableKinematicViscosity_)
{
insideCoeff_k += insideVolVars.kinematicViscosity();
outsideCoeff_k += outsideVolVars.kinematicViscosity();
insideCoeff_e += insideVolVars.kinematicViscosity();
outsideCoeff_e += outsideVolVars.kinematicViscosity();
}
Scalar insideCoeff_k = (insideVolVars.dynamicEddyViscosity() / insideVolVars.sigmaK()) + insideVolVars.viscosity();
Scalar outsideCoeff_k = (outsideVolVars.dynamicEddyViscosity() / outsideVolVars.sigmaK()) + outsideVolVars.viscosity();
Scalar insideCoeff_e = (insideVolVars.dynamicEddyViscosity() / insideVolVars.sigmaEpsilon()) + insideVolVars.viscosity();
Scalar outsideCoeff_e = (outsideVolVars.dynamicEddyViscosity() / outsideVolVars.sigmaEpsilon()) + outsideVolVars.viscosity();
// scale by extrusion factor
insideCoeff_k *= insideVolVars.extrusionFactor();
......
dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh
+ 6
− 6
View file @ 2a4fedcf
......@@ -89,8 +89,8 @@ public:
{
CellCenterPrimaryVariables storage = ParentType::computeStorageForCellCenter(problem, scv, volVars);
storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy();
storage[dissipationEqIdx] = volVars.dissipation();
storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy() * volVars.density();
storage[dissipationEqIdx] = volVars.dissipation() * volVars.density();
return storage;
}
......@@ -113,23 +113,23 @@ public:
// turbulence production is equal to dissipation -> exclude both terms (according to local equilibrium hypothesis, see FLUENT)
if (!volVars.isMatchingPoint())
{
source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.kinematicEddyViscosity()
source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.dynamicEddyViscosity()
* volVars.stressTensorScalarProduct();
}
source[dissipationEqIdx] += volVars.cOneEpsilon()
* dissipation / turbulentKineticEnergy
* 2.0 * volVars.kinematicEddyViscosity()
* 2.0 * volVars.dynamicEddyViscosity()
* volVars.stressTensorScalarProduct();
// destruction
// turbulence production is equal to dissipation -> exclude both terms (according to local equilibrium hypothesis, see FLUENT)
if (!volVars.isMatchingPoint())
{
source[turbulentKineticEnergyEqIdx] -= dissipation;
source[turbulentKineticEnergyEqIdx] -= dissipation * volVars.density();
}
source[dissipationEqIdx] -= volVars.cTwoEpsilon()
* dissipation * dissipation
/ turbulentKineticEnergy;
/ turbulentKineticEnergy * volVars.density();
return source;
}
......
dumux/freeflow/rans/twoeq/komega/model.hh
+ 12
− 12
View file @ 2a4fedcf
......@@ -30,33 +30,33 @@
*
* Turbulent Kinetic Energy balance:
* \f[
* \frac{\partial \left( k \right)}{\partial t}
* + \nabla \cdot \left( \mathbf{v} k \right)
* - \nabla \cdot \left[ \left( \nu + \sigma_\textrm{k} \nu_\textrm{t} \right) \nabla k \right]
* \frac{\partial \left( \varrho k \right)}{\partial t}
* + \nabla \cdot \left( \mathbf{v} \varrho k \right)
* - \nabla \cdot \left[ \left( \mu + \sigma_\textrm{k} \mu_\textrm{t} \right) \nabla k \right]
* - P
* + \beta_k^{*} k \omega
* + \beta_k^{*} k \varrho \omega
* = 0
* \f]
* with \f$ P = 2 \nu_\textrm{t} \mathbf{S} \cdot \mathbf{S} \f$
* with \f$ P = 2 \mu_\textrm{t} \mathbf{S} \cdot \mathbf{S} \f$
* and \f$ S_{ij} = \frac{1}{2} \left[ \frac{\partial}{\partial x_i} v_j + \frac{\partial}{\partial x_j} v_i \right] \f$
* based on \f$ a_{ij} \cdot b_{ij} = \sum_{i,j} a_{ij} b_{ij} \f$.
*
* Dissipation balance:
* \f[
* \frac{\partial \left( \omega \right)}{\partial t}
* + \nabla \cdot \left( \mathbf{v} \omega \right)
* - \nabla \cdot \left[ \left( \nu + \sigma_{\omega} \nu_\textrm{t} \right) \nabla \omega \right]
* \frac{\partial \left( \varrho \omega \right)}{\partial t}
* + \nabla \cdot \left( \mathbf{v} \varrho \omega \right)
* - \nabla \cdot \left[ \left( \mu + \sigma_{\omega} \mu_\textrm{t} \right) \nabla \omega \right]
* - \alpha \frac{\omega}{k} P
* + \beta_{\omega} \omega^2
* - \frac{\sigma_d}{\omega} \nabla k \nabla \omega
* - \varrho \frac{\sigma_d}{\omega} \nabla k \nabla \omega
* = 0
* \f]
*
* The kinematic eddy viscosity \f$ \nu_\textrm{t} \f$ is calculated as follows:
* \f[ \nu_\textrm{t} = \frac{k}{\tilde{\omega}} \f]
* The dynamic eddy viscosity \f$ \mu_\textrm{t} \f$ is calculated as follows:
* \f[ \mu_\textrm{t} = \varrho \frac{k}{\tilde{\omega}} \f]
*
* With a limited dissipation:
* \f[ \tilde{\omega} = \textrm{max} \left\{ \omega, 0.875 \sqrt{\frac{P}{\nu_\textrm{t} \beta_\textrm{k}}} \right\} \f]
* \f[ \tilde{\omega} = \textrm{max} \left\{ \omega, 0.875 \sqrt{\frac{P}{\mu_\textrm{t} \beta_\textrm{k}}} \right\} \f]
*
* And a cross-diffusion coefficient \f$ \sigma_\textrm{d} \f$
* \f[
......
dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh
+ 10
− 10
View file @ 2a4fedcf
......@@ -97,11 +97,11 @@ public:
// calculate advective flux
auto upwindTermK = [](const auto& volVars)
{
return volVars.turbulentKineticEnergy();
return volVars.turbulentKineticEnergy()*volVars.density();
};
auto upwindTermOmega = [](const auto& volVars)
{
return volVars.dissipation();
return volVars.dissipation()*volVars.density();
};
flux[turbulentKineticEnergyEqIdx]
......@@ -116,14 +116,14 @@ public:
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
// effective diffusion coefficients
Scalar insideCoeff_k = insideVolVars.kinematicViscosity()
+ ( insideVolVars.sigmaK() * insideVolVars.turbulentKineticEnergy() / insideVolVars.dissipation() );
Scalar outsideCoeff_k = outsideVolVars.kinematicViscosity()
+ ( outsideVolVars.sigmaK() * outsideVolVars.turbulentKineticEnergy() / outsideVolVars.dissipation() );
Scalar insideCoeff_w = insideVolVars.kinematicViscosity()
+ ( insideVolVars.sigmaOmega() * insideVolVars.turbulentKineticEnergy() / insideVolVars.dissipation() );
Scalar outsideCoeff_w = outsideVolVars.kinematicViscosity()
+ ( outsideVolVars.sigmaOmega() * outsideVolVars.turbulentKineticEnergy() / outsideVolVars.dissipation() );
Scalar insideCoeff_k = insideVolVars.viscosity()
+ ( insideVolVars.sigmaK() * insideVolVars.density()* insideVolVars.turbulentKineticEnergy() / insideVolVars.dissipation() );
Scalar outsideCoeff_k = outsideVolVars.viscosity()
+ ( outsideVolVars.sigmaK() * outsideVolVars.density()* outsideVolVars.turbulentKineticEnergy() / outsideVolVars.dissipation() );
Scalar insideCoeff_w = insideVolVars.viscosity()
+ ( insideVolVars.sigmaOmega() * insideVolVars.density()* insideVolVars.turbulentKineticEnergy() / insideVolVars.dissipation() );
Scalar outsideCoeff_w = outsideVolVars.viscosity()
+ ( outsideVolVars.sigmaOmega() * outsideVolVars.density()* outsideVolVars.turbulentKineticEnergy() / outsideVolVars.dissipation() );
// scale by extrusion factor
insideCoeff_k *= insideVolVars.extrusionFactor();
......
dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh
+ 7
− 7
View file @ 2a4fedcf
......@@ -83,8 +83,8 @@ public:
{
CellCenterPrimaryVariables storage = ParentType::computeStorageForCellCenter(problem, scv, volVars);
storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy();
storage[dissipationEqIdx] = volVars.dissipation();
storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy()*volVars.density();
storage[dissipationEqIdx] = volVars.dissipation()*volVars.density();
return storage;
}
......@@ -105,18 +105,18 @@ public:
// production
static const auto enableKOmegaProductionLimiter
= getParamFromGroup<bool>(problem.paramGroup(), "KOmega.EnableProductionLimiter", false);
Scalar productionTerm = 2.0 * volVars.kinematicEddyViscosity() * volVars.stressTensorScalarProduct();
Scalar productionTerm = 2.0 * volVars.dynamicEddyViscosity() * volVars.stressTensorScalarProduct();
if (enableKOmegaProductionLimiter)
{
Scalar productionAlternative = 20.0 * volVars.betaK() * volVars.turbulentKineticEnergy() * volVars.dissipation();
Scalar productionAlternative = 20.0 * volVars.density() * volVars.betaK() * volVars.turbulentKineticEnergy() * volVars.dissipation();
productionTerm = min(productionTerm, productionAlternative);
}
source[turbulentKineticEnergyEqIdx] += productionTerm;
source[dissipationEqIdx] += volVars.alpha() * volVars.dissipation() / volVars.turbulentKineticEnergy() * productionTerm;
// destruction
source[turbulentKineticEnergyEqIdx] -= volVars.betaK() * volVars.turbulentKineticEnergy() * volVars.dissipation();
source[dissipationEqIdx] -= volVars.betaOmega() * volVars.dissipation() * volVars.dissipation();
source[turbulentKineticEnergyEqIdx] -= volVars.betaK() * volVars.density() * volVars.turbulentKineticEnergy() * volVars.dissipation();
source[dissipationEqIdx] -= volVars.betaOmega() * volVars.density() * volVars.dissipation() * volVars.dissipation();
// cross-diffusion term
Scalar gradientProduct = 0.0;
......@@ -124,7 +124,7 @@ public:
gradientProduct += volVars.storedTurbulentKineticEnergyGradient()[i]
* volVars.storedDissipationGradient()[i];
if (gradientProduct > 0.0)
source[dissipationEqIdx] += 0.125 / volVars.dissipation() * gradientProduct;
source[dissipationEqIdx] += 0.125 *volVars.density() / volVars.dissipation() * gradientProduct;
return source;
}
......
dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh
+ 13
− 13
View file @ 2a4fedcf
......@@ -34,30 +34,30 @@
* \f$ \varepsilon = \tilde{\varepsilon} + D_\varepsilon \f$:
*
* \f[
* \frac{\partial \left( k \right)}{\partial t}
* + \nabla \cdot \left( \textbf{v} k \right)
* - \nabla \cdot \left( \left( \nu + \frac{\nu_\text{t}}{\sigma_\text{k}} \right) \nabla k \right)
* - 2 \nu_\text{t} \textbf{S} \cdot \textbf{S}
* + \tilde{\varepsilon}
* + D_\varepsilon
* \frac{\partial \left( \varrho k \right)}{\partial t}
* + \nabla \cdot \left( \textbf{v} \varhho k \right)
* - \nabla \cdot \left( \left( \mu + \frac{\mu_\text{t}}{\sigma_\text{k}} \right) \nabla k \right)
* - 2 \mu_\text{t} \textbf{S} \cdot \textbf{S}
* + \varrho \tilde{\varepsilon}
* + D_\varepsilon \varrho
* = 0
* \f].
*
* The dissipation balance is changed by introducing additional functions
* (\f$ E_\text{k}\f$, \f$ f_1 \f$, and \f$ f_2 \f$) to account for a dampening towards the wall:
* \f[
* \frac{\partial \left( \tilde{\varepsilon} \right)}{\partial t}
* + \nabla \cdot \left( \textbf{v} \tilde{\varepsilon} \right)
* - \nabla \cdot \left( \left( \nu + \frac{\nu_\text{t}}{\sigma_{\varepsilon}} \right) \nabla \tilde{\varepsilon} \right)
* - C_{1\tilde{\varepsilon}} f_1 \frac{\tilde{\varepsilon}}{k} 2 \nu_\text{t} \textbf{S} \cdot \textbf{S}
* + C_{2\tilde{\varepsilon}} f_2 \frac{\tilde{\varepsilon}^2}{k}
* - E_\text{k}
* \frac{\partial \left( \varrho \tilde{\varepsilon} \right)}{\partial t}
* + \nabla \cdot \left( \textbf{v} \varrho \tilde{\varepsilon} \right)
* - \nabla \cdot \left( \left( \mu + \frac{\mu_\text{t}}{\sigma_{\varepsilon}} \right) \nabla \tilde{\varepsilon} \right)
* - C_{1\tilde{\varepsilon}} f_1 \frac{\tilde{\varepsilon}}{k} 2 \mu_\text{t} \textbf{S} \cdot \textbf{S}
* + C_{2\tilde{\varepsilon}} \varrho f_2 \frac{\tilde{\varepsilon}^2}{k}
* - E_\text{k} \varrho
* = 0
* \f].
*
* The kinematic eddy viscosity \f$ \nu_\text{t} \f$ is dampened by \f$ f_\mu \f$:
* \f[
* \nu_\text{t} = C_\mu f_\mu \frac{k^2}{\tilde{\varepsilon}}
* \mu_\text{t} = \varrho C_\mu f_\mu \frac{k^2}{\tilde{\varepsilon}}
* \f].
*
* The auxiliary and dampening functions are defined as:
......
dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
+ 10
− 10
View file @ 2a4fedcf
......@@ -97,11 +97,11 @@ public:
// calculate advective flux
auto upwindTermK = [](const auto& volVars)
{
return volVars.turbulentKineticEnergy();
return volVars.turbulentKineticEnergy() * volVars.density();
};
auto upwindTermEpsilon = [](const auto& volVars)
{
return volVars.dissipationTilde();
return volVars.dissipationTilde() * volVars.density();
};
flux[turbulentKineticEnergyEqIdx]
......@@ -116,14 +116,14 @@ public:
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
// effective diffusion coefficients
Scalar insideCoeff_k = insideVolVars.kinematicViscosity()
+ insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaK();
Scalar outsideCoeff_k = outsideVolVars.kinematicViscosity()
+ outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaK();
Scalar insideCoeff_e = insideVolVars.kinematicViscosity()
+ insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaEpsilon();
Scalar outsideCoeff_e = outsideVolVars.kinematicViscosity()
+ outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaEpsilon();
Scalar insideCoeff_k = insideVolVars.viscosity() + insideVolVars.kinematicEddyViscosity()
* insideVolVars.density() / insideVolVars.sigmaK();
Scalar outsideCoeff_k = outsideVolVars.viscosity() + outsideVolVars.kinematicEddyViscosity()
* outsideVolVars.density() / outsideVolVars.sigmaK();
Scalar insideCoeff_e = insideVolVars.viscosity() + insideVolVars.kinematicEddyViscosity()
* insideVolVars.density() / insideVolVars.sigmaEpsilon();
Scalar outsideCoeff_e = outsideVolVars.viscosity() + outsideVolVars.kinematicEddyViscosity()
* outsideVolVars.density() / outsideVolVars.sigmaEpsilon();
// scale by extrusion factor
insideCoeff_k *= insideVolVars.extrusionFactor();
......
dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/localresidual.hh
+ 8
− 8
View file @ 2a4fedcf
......@@ -82,8 +82,8 @@ public:
{
CellCenterPrimaryVariables storage = ParentType::computeStorageForCellCenter(problem, scv, volVars);
storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy();
storage[dissipationEqIdx] = volVars.dissipationTilde();
storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy() * volVars.density();
storage[dissipationEqIdx] = volVars.dissipationTilde() * volVars.density();
return storage;
}
......@@ -101,22 +101,22 @@ public:
const auto& volVars = elemVolVars[scv];
// production
source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.kinematicEddyViscosity()
source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.dynamicEddyViscosity()
* volVars.stressTensorScalarProduct();
source[dissipationEqIdx] += volVars.cOneEpsilon() * volVars.fOne()
* volVars.dissipationTilde() / volVars.turbulentKineticEnergy()
* 2.0 * volVars.kinematicEddyViscosity()
* 2.0 * volVars.dynamicEddyViscosity()
* volVars.stressTensorScalarProduct();
// destruction
source[turbulentKineticEnergyEqIdx] -= volVars.dissipationTilde();
source[dissipationEqIdx] -= volVars.cTwoEpsilon() * volVars.fTwo()
source[turbulentKineticEnergyEqIdx] -= volVars.dissipationTilde() * volVars.density();
source[dissipationEqIdx] -= volVars.cTwoEpsilon() * volVars.fTwo() * volVars.density()
* volVars.dissipationTilde() * volVars.dissipationTilde()
/ volVars.turbulentKineticEnergy();
// dampening functions
source[turbulentKineticEnergyEqIdx] -= volVars.dValue();
source[dissipationEqIdx] += volVars.eValue();
source[turbulentKineticEnergyEqIdx] -= volVars.dValue() * volVars.density();
source[dissipationEqIdx] += volVars.eValue() * volVars.density();
return source;
}
......
test/freeflow/rans/main.cc
+ 1
− 1
View file @ 2a4fedcf
......@@ -45,7 +45,7 @@
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/discretization/method.hh>
#include "problem.hh"
#include "properties.hh"
/*!
* \brief Provides an interface for customizing error messages associated with
......
test/freeflow/rans/problem.hh
+ 2
− 75
View file @ 2a4fedcf
......@@ -27,85 +27,12 @@
#ifndef DUMUX_PIPE_LAUFER_PROBLEM_HH
#define DUMUX_PIPE_LAUFER_PROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
#include <dumux/freeflow/turbulenceproperties.hh>
#include <dumux/freeflow/turbulencemodel.hh>
#include <dumux/common/properties.hh>
#include <dumux/freeflow/rans/problem.hh>
#include <dumux/freeflow/rans/oneeq/problem.hh>
#include <dumux/freeflow/rans/oneeq/model.hh>
#include <dumux/freeflow/rans/twoeq/kepsilon/model.hh>
#include <dumux/freeflow/rans/twoeq/kepsilon/problem.hh>
#include <dumux/freeflow/rans/twoeq/komega/model.hh>
#include <dumux/freeflow/rans/twoeq/komega/problem.hh>
#include <dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh>
#include <dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh>
#include <dumux/freeflow/rans/zeroeq/model.hh>
#include <dumux/material/components/air.hh>
#include <dumux/material/fluidsystems/1pgas.hh>
#include <dumux/freeflow/turbulencemodel.hh>
namespace Dumux {
template <class TypeTag>
class PipeLauferProblem;
namespace Properties {
// Create new type tags
namespace TTag {
// Base Typetag
struct RANSModel { using InheritsFrom = std::tuple<StaggeredFreeFlowModel>; };
// Isothermal Typetags
struct PipeLauferZeroEq { using InheritsFrom = std::tuple<RANSModel, ZeroEq>; };
struct PipeLauferOneEq { using InheritsFrom = std::tuple<RANSModel, OneEq>; };
struct PipeLauferKOmega { using InheritsFrom = std::tuple<RANSModel, KOmega>; };
struct PipeLauferLowReKEpsilon { using InheritsFrom = std::tuple<RANSModel, LowReKEpsilon>; };
struct PipeLauferKEpsilon { using InheritsFrom = std::tuple<RANSModel, KEpsilon>; };
// Non-Isothermal Typetags
struct PipeLauferNIZeroEq { using InheritsFrom = std::tuple<RANSModel, ZeroEqNI>; };
struct PipeLauferNIOneEq { using InheritsFrom = std::tuple<RANSModel, OneEqNI>; };
struct PipeLauferNIKOmega { using InheritsFrom = std::tuple<RANSModel, KOmegaNI>; };
struct PipeLauferNILowReKEpsilon { using InheritsFrom = std::tuple<RANSModel, LowReKEpsilonNI>; };
struct PipeLauferNIKEpsilon { using InheritsFrom = std::tuple<RANSModel, KEpsilonNI>; };
} // end namespace TTag
// the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::RANSModel>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = FluidSystems::OnePGas<Scalar, Components::Air<Scalar> >;
};
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::RANSModel>
{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::RANSModel>
{ using type = Dumux::PipeLauferProblem<TypeTag>; };
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::RANSModel> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::RANSModel> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::RANSModel> { static constexpr bool value = true; };
} // end namespace Properties
/*!
* \ingroup NavierStokesTests
* \brief Test problem for the one-phase (Navier-) Stokes problem in a channel.
......
test/freeflow/rans/properties.hh 0 → 100644
+ 98
− 0
View file @ 2a4fedcf
// -*- 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 3 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
* \ingroup RANSTests
* \brief Pipe flow test for the staggered grid RANS model
*
* This test simulates pipe flow experiments performed by John Laufer in 1954
* \cite Laufer1954a.
*/
#ifndef DUMUX_PIPE_LAUFER_PROPERTIES_HH
#define DUMUX_PIPE_LAUFER_PROPERTIES_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/material/components/air.hh>
#include <dumux/material/fluidsystems/1pgas.hh>
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
#include <dumux/freeflow/turbulenceproperties.hh>
#include <dumux/freeflow/turbulencemodel.hh>
#include <dumux/freeflow/rans/problem.hh>
#include <dumux/freeflow/rans/zeroeq/model.hh>
#include <dumux/freeflow/rans/oneeq/problem.hh>
#include <dumux/freeflow/rans/oneeq/model.hh>
#include <dumux/freeflow/rans/twoeq/kepsilon/model.hh>
#include <dumux/freeflow/rans/twoeq/kepsilon/problem.hh>
#include <dumux/freeflow/rans/twoeq/komega/model.hh>
#include <dumux/freeflow/rans/twoeq/komega/problem.hh>
#include <dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh>
#include <dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh>
#include "problem.hh"
namespace Dumux::Properties {
// Create new type tags
namespace TTag {
// Base Typetag
struct RANSModel { using InheritsFrom = std::tuple<StaggeredFreeFlowModel>; };
// Isothermal Typetags
struct PipeLauferZeroEq { using InheritsFrom = std::tuple<RANSModel, ZeroEq>; };
struct PipeLauferOneEq { using InheritsFrom = std::tuple<RANSModel, OneEq>; };
struct PipeLauferKOmega { using InheritsFrom = std::tuple<RANSModel, KOmega>; };
struct PipeLauferLowReKEpsilon { using InheritsFrom = std::tuple<RANSModel, LowReKEpsilon>; };
struct PipeLauferKEpsilon { using InheritsFrom = std::tuple<RANSModel, KEpsilon>; };
// Non-Isothermal Typetags
struct PipeLauferNIZeroEq { using InheritsFrom = std::tuple<RANSModel, ZeroEqNI>; };
struct PipeLauferNIOneEq { using InheritsFrom = std::tuple<RANSModel, OneEqNI>; };
struct PipeLauferNIKOmega { using InheritsFrom = std::tuple<RANSModel, KOmegaNI>; };
struct PipeLauferNILowReKEpsilon { using InheritsFrom = std::tuple<RANSModel, LowReKEpsilonNI>; };
struct PipeLauferNIKEpsilon { using InheritsFrom = std::tuple<RANSModel, KEpsilonNI>; };
} // end namespace TTag
// the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::RANSModel>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = FluidSystems::OnePGas<Scalar, Components::Air<Scalar> >;
};
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::RANSModel>
{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::RANSModel>
{ using type = Dumux::PipeLauferProblem<TypeTag>; };
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::RANSModel> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::RANSModel> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::RANSModel> { static constexpr bool value = true; };
} // end namespace Dumux::Properties
#endif
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