From a444e7b35e5b2c7db571643bb2eb1e94d3007d2b Mon Sep 17 00:00:00 2001
From: Thomas Fetzer <thomas.fetzer@iws.uni-stuttgart.de>
Date: Thu, 13 Feb 2014 13:12:15 +0000
Subject: [PATCH] [stokes] implement switch for symmetrizing of the viscous
term, symmetrizing is default FS#212. Deprecated eddyViscosity().
reviewed by klaus
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@12468 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/freeflow/stokes/stokesfluxvariables.hh | 18 +++-
dumux/freeflow/stokes/stokeslocalresidual.hh | 98 +++++++++++++++----
dumux/freeflow/stokes/stokesproperties.hh | 1 +
.../freeflow/stokes/stokespropertydefaults.hh | 5 +-
test/freeflow/stokes/stokestestproblem.hh | 5 +-
5 files changed, 103 insertions(+), 24 deletions(-)
diff --git a/dumux/freeflow/stokes/stokesfluxvariables.hh b/dumux/freeflow/stokes/stokesfluxvariables.hh
index 3dd49e83b8..ba5925fb84 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 c60bd5c7f4..aa8843562a 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 f6139ea91a..8298815e42 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 0ecaddb678..ab64eb0999 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 33016829c3..02a1409230 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);
}
--
GitLab