Commit a444e7b3 authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[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
parent 5de8c537
......@@ -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.
......
......@@ -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);
}
}
}
......
......@@ -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
......
......@@ -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);
}
......
......@@ -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);
}
......
Supports Markdown
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