diff --git a/dumux/freeflow/stokes/stokesfluxvariables.hh b/dumux/freeflow/stokes/stokesfluxvariables.hh
index 3dd49e83b842b10a2bd95662adc0740a38bde9dc..ba5925fb84dd16f61b1f3dd8d3847f23ff23653c 100644
--- a/dumux/freeflow/stokes/stokesfluxvariables.hh
+++ b/dumux/freeflow/stokes/stokesfluxvariables.hh
@@ -231,11 +231,27 @@ public:
/*!
* \brief Return the dynamic eddy viscosity
- * \f$\mathrm{[Pa\cdot s]} = \mathrm{[N\cdot s/m^2]}\f$ (if implemented).
+ * \f$\mathrm{[Pa \cdot s]} = \mathrm{[N \cdot s/m^2]}\f$ (if implemented).
*/
+ DUNE_DEPRECATED_MSG("Function eddyViscosity() is deprecated, use dynamicEddyViscosity() instead.")
const Scalar eddyViscosity() const
+ { return dynamicEddyViscosity(); }
+
+ /*!
+ * \brief Return the dynamic eddy viscosity
+ * \f$\mathrm{[Pa \cdot s]} = \mathrm{[N \cdot s/m^2]}\f$ (if implemented).
+ */
+ const Scalar dynamicEddyViscosity() const
+ { return kinematicEddyViscosity() * density(); }
+
+ /*!
+ * \brief Return the kinematic eddy viscosity
+ * \f$\mathrm{[m^2/s]}\f$ (if implemented).
+ */
+ const Scalar kinematicEddyViscosity() const
{ return 0; }
+
// /*!
// * \brief Return the divergence of the normal velocity at the
// * integration point.
diff --git a/dumux/freeflow/stokes/stokeslocalresidual.hh b/dumux/freeflow/stokes/stokeslocalresidual.hh
index c60bd5c7f408643a7c70996395061e84df1dee3a..aa8843562aa313ad0c4bd2dcaa2716d473bf2b34 100644
--- a/dumux/freeflow/stokes/stokeslocalresidual.hh
+++ b/dumux/freeflow/stokes/stokeslocalresidual.hh
@@ -90,7 +90,8 @@ protected:
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
-
+
+ static const bool enableUnsymmetrizedVelocityGradient = GET_PROP_VALUE(TypeTag, EnableUnsymmetrizedVelocityGradient);
static const bool calculateNavierStokes = GET_PROP_VALUE(TypeTag, EnableNavierStokes);
public:
@@ -215,22 +216,30 @@ protected:
// at the center of the SCV in computeSource
// dynamic viscosity is upwinded
- // compute symmetrized gradient for the momentum flux:
- // mu (grad v + (grad v)^t)
- Dune::FieldMatrix<Scalar, dim, dim> symmVelGrad = fluxVars.velocityGrad();
- for (int i=0; i<dim; ++i)
- for (int j=0; j<dim; ++j)
- symmVelGrad[i][j] += fluxVars.velocityGrad()[j][i];
+ Dune::FieldMatrix<Scalar, dim, dim> velGrad = fluxVars.velocityGrad();
+ if (enableUnsymmetrizedVelocityGradient)
+ {
+ // nothing has to be done in this case:
+ // grad v
+ }
+ else
+ {
+ // compute symmetrized gradient for the momentum flux:
+ // grad v + (grad v)^T
+ for (int i=0; i<dim; ++i)
+ for (int j=0; j<dim; ++j)
+ velGrad[i][j] += fluxVars.velocityGrad()[j][i];
+ }
DimVector velGradComp(0.);
for (int velIdx = 0; velIdx < dim; ++velIdx)
{
- velGradComp = symmVelGrad[velIdx];
+ velGradComp = velGrad[velIdx];
// TODO: dilatation term has to be accounted for in outflow, coupling, neumann
// velGradComp[velIdx] += 2./3*fluxVars.velocityDiv;
- velGradComp *= fluxVars.dynamicViscosity() + fluxVars.eddyViscosity();
+ velGradComp *= fluxVars.dynamicViscosity() + fluxVars.dynamicEddyViscosity();
flux[momentumXIdx + velIdx] -=
velGradComp*fluxVars.face().normal;
@@ -244,7 +253,7 @@ protected:
// gravityTerm;
}
-
+
// this term changes the Stokes equation to the Navier-Stokes equation
// rho v (v*n)
// rho and first v are upwinded, second v is evaluated at the face
@@ -380,12 +389,27 @@ protected:
// the fluxes at the boundary are added in the second step
if (momentumBalanceDirichlet_(bcTypes))
{
- DimVector muGradVelNormal(0.);
const DimVector &boundaryFaceNormal =
boundaryVars.face().normal;
- boundaryVars.velocityGrad().umv(boundaryFaceNormal, muGradVelNormal);
- muGradVelNormal *= boundaryVars.dynamicViscosity();
+ Dune::FieldMatrix<Scalar, dim, dim> velGrad = boundaryVars.velocityGrad();
+ if (enableUnsymmetrizedVelocityGradient)
+ {
+ // nothing has to be done in this case:
+ // grad v
+ }
+// else
+// {
+// // compute symmetrized gradient for the momentum flux:
+// // grad v + (grad v)^T
+// for (int i=0; i<dim; ++i)
+// for (int j=0; j<dim; ++j)
+// velGrad[i][j] += boundaryVars.velocityGrad()[j][i];
+// }
+
+ DimVector muGradVelNormal(0.0);
+ velGrad.umv(boundaryFaceNormal, muGradVelNormal);
+ muGradVelNormal *= boundaryVars.dynamicViscosity() + boundaryVars.dynamicEddyViscosity();
for (unsigned int i=0; i < this->residual_.size(); i++)
Valgrind::CheckDefined(this->residual_[i]);
@@ -467,14 +491,30 @@ protected:
const DimVector& elementUnitNormal = isIt->centerUnitOuterNormal();
tangent[0] = elementUnitNormal[1]; //TODO: 3D
tangent[1] = -elementUnitNormal[0];
- DimVector tangentialVelGrad;
- boundaryVars.velocityGrad().mv(tangent, tangentialVelGrad);
- tangentialVelGrad *= boundaryVars.dynamicViscosity();
+
+ Dune::FieldMatrix<Scalar, dim, dim> velGrad = boundaryVars.velocityGrad();
+ if (enableUnsymmetrizedVelocityGradient)
+ {
+ // nothing has to be done in this case:
+ // grad v
+ }
+// else
+// {
+// // compute symmetrized gradient for the momentum flux:
+// // mu (grad v + (grad v)^t)
+// for (int i=0; i<dim; ++i)
+// for (int j=0; j<dim; ++j)
+// velGrad[i][j] += boundaryVars.velocityGrad()[j][i];
+// }
+
+ DimVector muGradVelTangential(0.0);
+ velGrad.mv(tangent, muGradVelTangential);
+ muGradVelTangential *= boundaryVars.dynamicViscosity() + boundaryVars.dynamicEddyViscosity();
this->residual_[scvIdx][massBalanceIdx] -= stabilizationBeta_*0.5*
this->curVolVars_(scvIdx).pressure();
this->residual_[scvIdx][massBalanceIdx] -= stabilizationBeta_*0.5*
- (tangentialVelGrad*tangent);
+ (muGradVelTangential*tangent);
for (int momentumIdx = momentumXIdx; momentumIdx <= lastMomentumIdx; momentumIdx++)
this->residual_[scvIdx][massBalanceIdx] -= stabilizationBeta_*0.5
@@ -525,12 +565,28 @@ protected:
const DimVector& elementUnitNormal = isIt->centerUnitOuterNormal();
tangent[0] = elementUnitNormal[1];
tangent[1] = -elementUnitNormal[0];
- DimVector tangentialVelGrad;
- boundaryVars.velocityGrad().mv(tangent, tangentialVelGrad);
+ Dune::FieldMatrix<Scalar, dim, dim> velGrad = boundaryVars.velocityGrad();
+ if (enableUnsymmetrizedVelocityGradient)
+ {
+ // nothing has to be done in this case:
+ // grad v
+ }
+ else
+ {
+ // compute symmetrized gradient for the momentum flux:
+ // mu (grad v + (grad v)^t)
+ for (int i=0; i<dim; ++i)
+ for (int j=0; j<dim; ++j)
+ velGrad[i][j] += boundaryVars.velocityGrad()[j][i];
+ }
+
+ DimVector muGradVelTangential(0.0);
+ velGrad.umv(tangent, muGradVelTangential);
+ muGradVelTangential *= boundaryVars.dynamicViscosity() + boundaryVars.dynamicEddyViscosity();
+
this->residual_[scvIdx][massBalanceIdx] -= 0.5*stabilizationBeta_
- * boundaryVars.dynamicViscosity()
- * (tangentialVelGrad*tangent);
+ * (muGradVelTangential*tangent);
}
}
}
diff --git a/dumux/freeflow/stokes/stokesproperties.hh b/dumux/freeflow/stokes/stokesproperties.hh
index f6139ea91a71bf4829fafc680f5cc2bc52bc8871..8298815e427702f105b3f51166848b7042f89bff 100644
--- a/dumux/freeflow/stokes/stokesproperties.hh
+++ b/dumux/freeflow/stokes/stokesproperties.hh
@@ -55,6 +55,7 @@ NEW_PROP_TAG(FluidSystem); //!< The employed fluid system
NEW_PROP_TAG(FluidState);
NEW_PROP_TAG(StokesStabilizationAlpha); //!< The parameter for the stabilization
NEW_PROP_TAG(StokesStabilizationBeta); //!< The parameter for the stabilization at boundaries
+NEW_PROP_TAG(EnableUnsymmetrizedVelocityGradient); //!< Returns whether unsymmetrized velocity gradient for viscous term is used
NEW_PROP_TAG(EnableNavierStokes); //!< Returns whether Navier-Stokes should be solved instead of plain Stokes
NEW_PROP_TAG(PhaseIdx); //!< A phase index in case that a two-phase fluidsystem is used
diff --git a/dumux/freeflow/stokes/stokespropertydefaults.hh b/dumux/freeflow/stokes/stokespropertydefaults.hh
index 0ecaddb678f2b5f8588f150553b4d3a4620a0817..ab64eb09999f4fccb1e2ba647a9676fa9f664cf0 100644
--- a/dumux/freeflow/stokes/stokespropertydefaults.hh
+++ b/dumux/freeflow/stokes/stokespropertydefaults.hh
@@ -134,6 +134,9 @@ SET_BOOL_PROP(BoxStokes, EvalGradientsAtSCVCenter, true);
//! Set the phaseIndex per default to zero (important for two-phase fluidsystems).
SET_INT_PROP(BoxStokes, PhaseIdx, 0);
+//! Use symmetrizedVelocityGradient by default
+SET_BOOL_PROP(BoxStokes, EnableUnsymmetrizedVelocityGradient, false);
+
//! Set calculation to Stokes, not Navier-Stokes
SET_BOOL_PROP(BoxStokes, EnableNavierStokes, false);
@@ -143,7 +146,7 @@ SET_SCALAR_PROP(BoxStokes, StokesStabilizationAlpha, 0.0);
//! Stabilization factor for the boundaries
SET_SCALAR_PROP(BoxStokes, StokesStabilizationBeta, 0.0);
-// enable gravity by default
+//! Set gravity by default
SET_BOOL_PROP(BoxStokes, ProblemEnableGravity, true);
}
diff --git a/test/freeflow/stokes/stokestestproblem.hh b/test/freeflow/stokes/stokestestproblem.hh
index 33016829c317bb8a03ee7cdab0d6f08f5f087e0b..02a1409230efe3f9c4d4d06d83a175e2c9fa8276 100644
--- a/test/freeflow/stokes/stokestestproblem.hh
+++ b/test/freeflow/stokes/stokestestproblem.hh
@@ -67,7 +67,10 @@ SET_TYPE_PROP(BoxStokes, Scalar, double);
//! A stabilization factor. Set negative for stabilization and to zero for no stabilization
SET_SCALAR_PROP(StokesTestProblem, StokesStabilizationAlpha, -1.0);
-// Enable gravity
+//! Disable unsymmetrized velecoity gradient by default
+SET_BOOL_PROP(StokesTestProblem, EnableUnsymmetrizedVelocityGradient, false);
+
+//! Disable gravity by default
SET_BOOL_PROP(StokesTestProblem, ProblemEnableGravity, false);
}