Commit afef51b0 authored by Martin Utz's avatar Martin Utz
Browse files

[swe][frictionlaws] Return shear stress instead of source

The friction law compute now the bed shear stress instead of the source
term. So they return a quantity with a clear physical meaning, which can
be also used for other purpuses (e.g. sediment transport). Therefore the
computeSource method was renamed to computeShearStress.
computeShearStress uses now volvars instead of the explicit use of
waterDepth, u, v.
parent 241c5038
......@@ -35,24 +35,21 @@ namespace Dumux {
* The LET mobility model is used to limit the friction for small water depths.
*/
template <typename NumEqVector>
template <typename Scalar, typename VolumeVariables >
class FrictionLaw
{
using Scalar = typename NumEqVector::value_type;
public:
/*!
* \brief Compute the friction source term.
* \brief Compute the shear stress.
*
* The friction source term contains the losses due to friction.
* Compute the shear stress due to friction. The shear stress is not a tensor as know
* from contiuums mechanics, but a force projected on an area. Therefore it is a
* vector with two entries.
*
* \return Friction source term.
* \return shear stress. First entry is the x-component, the second the y-component.
*/
virtual NumEqVector computeSource(const Scalar waterDepth,
const Scalar frictionValue,
const Scalar u,
const Scalar v) const = 0;
virtual Dune::FieldVector<Scalar, 2> computeShearStress(const VolumeVariables& VolVar) const = 0;
/*!
* \brief Limit the friction for small water depth.
......
......@@ -36,48 +36,50 @@ namespace Dumux {
* The LET mobility model is used to limit the friction for small water depths.
*/
template <typename NumEqVector>
class FrictionLawManning : public FrictionLaw<NumEqVector>
template <typename Scalar, typename VolumeVariables>
class FrictionLawManning : public FrictionLaw<Scalar, VolumeVariables>
{
using Scalar = typename NumEqVector::value_type;
public:
FrictionLawManning(const Scalar gravity)
: gravity_(gravity) {}
/*!
* \brief Constructor
*
* \param gravity Gravity constant [m/s^2]
* \param manningN Manning friction coefficient [-]
*/
FrictionLawManning(const Scalar gravity, const Scalar manningN)
: gravity_(gravity), manningN_(manningN) {}
/*!
* \brief Compute the friction source term.
* \brief Compute the shear stress.
*
* \param volVars Volume variables
*
* \param waterDepth water depth.
* \param frictionValue Manning friction value.
* \param u velocity in x-direction.
* \param v velocity in y-direction.
* Compute the shear stress due to friction. The shear stress is not a tensor as know
* from contiuums mechanics, but a force projected on an area. Therefore it is a
* vector with two entries.
*
* \return Friction source term.
* \return shear stress. First entry is the x-component, the second the y-component.
*/
NumEqVector computeSource(const Scalar waterDepth,
const Scalar frictionValue,
const Scalar u,
const Scalar v) const final
Dune::FieldVector<Scalar, 2> computeShearStress(const VolumeVariables& volVars) const final
{
using std::pow;
using std::hypot;
NumEqVector source(0.0);
Dune::FieldVector<Scalar, 2> shearStress(0.0);
Scalar roughnessHeight = pow(25.68/(1.0/frictionValue),6.0);
roughnessHeight = this->limitRoughH(roughnessHeight, waterDepth);
const Scalar c = pow((waterDepth + roughnessHeight),1.0/6.0) * 1.0/(frictionValue);
const Scalar uv = hypot(u,v);
Scalar roughnessHeight = pow(25.68/(1.0/manningN_),6.0);
roughnessHeight = this->limitRoughH(roughnessHeight, volVars.waterDepth());
const Scalar c = pow((volVars.waterDepth() + roughnessHeight),1.0/6.0) * 1.0/(manningN_);
const Scalar uv = hypot(volVars.velocity(0),volVars.velocity(1));
source[0] = 0.0;
source[1] = -gravity_/(c*c) * u * uv;
source[2] = -gravity_/(c*c) * v * uv;
shearStress[0] = -gravity_/(c*c) * volVars.velocity(0) * uv;
shearStress[1] = -gravity_/(c*c) * volVars.velocity(1) * uv;
return source;
return shearStress;
}
private:
Scalar gravity_;
Scalar manningN_;
};
} // end namespace Dumux
......
......@@ -36,45 +36,49 @@ namespace Dumux {
* The LET mobility model is used to limit the friction for small water depths.
*/
template <typename NumEqVector>
class FrictionLawNikuradse : public FrictionLaw<NumEqVector>
template <typename Scalar, typename VolumeVariables>
class FrictionLawNikuradse : public FrictionLaw<Scalar, VolumeVariables>
{
using Scalar = typename NumEqVector::value_type;
public:
/*!
* \brief Compute the friction source term.
*
* \param waterDepth water depth.
* \param frictionValue The equivalent sand roughness.
* \param u velocity in x-direction.
* \param v velocity in y-direction.
* \brief Constructor
*
* \return Friction source term.
* \param ks Equivalent sand roughness [m]
*/
FrictionLawNikuradse(const Scalar ks)
: ks_(ks) {}
NumEqVector computeSource(const Scalar waterDepth,
const Scalar frictionValue,
const Scalar u,
const Scalar v) const final
/*!
* \brief Compute the shear stress.
*
* \param volVars Volume variables
*
* Compute the shear stress due to friction. The shear stress is not a tensor as know
* from contiuums mechanics, but a force projected on an area. Therefore it is a
* vector with two entries.
*
* \return shear stress. First entry is the x-component, the second the y-component.
*/
Dune::FieldVector<Scalar, 2> computeShearStress(const VolumeVariables& volVars) const final
{
using std::pow;
using std::log;
using std::hypot;
NumEqVector source(0.0);
Dune::FieldVector<Scalar, 2> shearStress(0.0);
Scalar roughnessHeight = frictionValue;
roughnessHeight = this->limitRoughH(roughnessHeight, waterDepth);
const Scalar ustarH = pow(0.41,2.0)/pow(log((12*(waterDepth + roughnessHeight))/frictionValue),2.0);
const Scalar uv = hypot(u,v);
Scalar roughnessHeight = ks_;
roughnessHeight = this->limitRoughH(roughnessHeight, volVars.waterDepth());
const Scalar ustarH = pow(0.41,2.0)/pow(log((12*(volVars.waterDepth() + roughnessHeight))/ks_),2.0);
const Scalar uv = hypot(volVars.velocity(0),volVars.velocity(1));
source[0] = 0.0;
source[1] = -ustarH * u * uv;
source[2] = -ustarH * v * uv;
shearStress[0] = -ustarH * volVars.velocity(0) * uv;
shearStress[1] = -ustarH * volVars.velocity(1) * uv;
return source;
return shearStress;
}
private:
Scalar ks_;
};
} // end namespace Dumux
......
......@@ -124,6 +124,7 @@ class RoughChannelProblem : public ShallowWaterProblem<TypeTag>
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using NeumannFluxes = GetPropType<TypeTag, Properties::NumEqVector>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = GetPropType<TypeTag, Properties::GridView>;
......@@ -147,11 +148,11 @@ public:
const auto gravity = this->spatialParams().gravity();
if (getParam<std::string>("Problem.FrictionLaw") == ("Manning"))
{
frictionLaw_ = std::make_shared<FrictionLawManning<NumEqVector>>(gravity);
frictionLaw_ = std::make_shared<FrictionLawManning<Scalar, VolumeVariables>>(gravity, constManningN_);
}
else if (getParam<std::string>("Problem.FrictionLaw") == ("Nikuradse"))
{
frictionLaw_ = std::make_shared<FrictionLawNikuradse<NumEqVector>>();
frictionLaw_ = std::make_shared<FrictionLawNikuradse<Scalar, VolumeVariables>>(0.003);
std::cout<<"\nWARNING: This test is meant to be run for the friction law after Manning.\n";
std::cout<<" You are running it with Nikuradse, although the friction values in\n";
std::cout<<" the input file and the analytic solution don't fit to Nikuradse!\n\n.";
......@@ -238,20 +239,34 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
using std::hypot;
const auto& globalPos = scv.center();
const auto& volVars = elemVolVars[scv];
const Scalar manningN = this->spatialParams().frictionValue(globalPos);
NumEqVector source (0.0);
const Scalar h = volVars.waterDepth();
const Scalar u = volVars.velocity(0);
const Scalar v = volVars.velocity(1);
NumEqVector source = frictionLaw_->computeSource(h, manningN, u, v);
Dune::FieldVector<Scalar, 2> bottomShearStress = frictionLaw_->computeShearStress(volVars);
source += computeBottomFrictionSource(bottomShearStress);
return source;
}
/*!
* \brief Compute the source term due to bottom friction
*
* \param bottomShearStress Shear stress due to bottom friction.
*
* \return source
*/
NumEqVector computeBottomFrictionSource(const Dune::FieldVector<Scalar, 2> bottomShearStress) const
{
NumEqVector bottomFrictionSource(0.0);
bottomFrictionSource[0] = 0.0;
bottomFrictionSource[1] =bottomShearStress[0];
bottomFrictionSource[2] =bottomShearStress[1];
return bottomFrictionSource;
}
// \}
......@@ -380,7 +395,7 @@ private:
Scalar constManningN_; // analytic solution is only available for const friction.
Scalar bedSlope_;
Scalar discharge_; // discharge at the inflow boundary
std::shared_ptr<FrictionLaw<NumEqVector>> frictionLaw_;
std::shared_ptr<FrictionLaw<Scalar, VolumeVariables>> frictionLaw_;
static constexpr Scalar eps_ = 1.0e-6;
std::string name_;
};
......
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