Commit 0884e520 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

Merge branch 'feature/free-flow-pseudo3d' into 'master'

[navierstokes] Add convenience function to calculate pseudo 3D friction term

See merge request !925
parents ec6f5f8e 5cd8a129
......@@ -1068,7 +1068,7 @@ url = {http://www.sciencedirect.com/science/article/pii/S0169772204001160}
volume = {198},
pages = {71--78}
}
@misc{cooper2008,
@misc{cooper2008,
title={{Release of the IAPWS formulation 2008 for the viscosity of ordinary water substance}},
author={Cooper, J. R. and Dooley, R. B.},
year={2008},
......@@ -1676,3 +1676,29 @@ url={http://dx.doi.org/10.1007/s11242-015-0599-1}
doi = {10.1115/1.3240815},
url = {http://dx.doi.org/10.1115/1.3240815},
}
@Article{flekkoy1995a,
title={{Hydrodynamic dispersion at stagnation points: Simulations and experiments}},
author={{Flekk{\o}y, EG and Oxaal, U and Feder, J and J{\o}ssang, T}},
journal={Physical Review E},
volume={52},
number={5},
pages={4952},
year={1995},
publisher={APS},
doi = {10.1103/PhysRevE.52.4952},
url = {https://doi.org/10.1103/PhysRevE.52.4952}
}
@Article{venturoli2006a,
title={{Two-dimensional lattice-Boltzmann simulations of single phase flow in a pseudo two-dimensional micromodel}},
author={{Venturoli, Maddalena and Boek, Edo S}},
journal={Physica A: Statistical Mechanics and its Applications},
volume={362},
number={1},
pages={23--29},
year={2006},
publisher={Elsevier},
doi = {10.1016/j.physa.2005.09.006},
url = {https://doi.org/10.1016/j.physa.2005.09.006}
}
......@@ -116,19 +116,57 @@ public:
* If the <tt>Problem.EnableGravity</tt> parameter is true, this means
* \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$
*/
const GlobalPosition &gravity() const
const GlobalPosition& gravity() const
{ return gravity_; }
//! Applys the initial face solution (velocities on the faces). Specialization for staggered grid discretization.
template <class G = FVGridGeometry>
typename std::enable_if<G::discMethod == DiscretizationMethod::staggered, void>::type
applyInitialFaceSolution(SolutionVector& sol,
const SubControlVolumeFace& scvf,
const PrimaryVariables& initSol) const
const SubControlVolumeFace& scvf,
const PrimaryVariables& initSol) const
{
sol[FVGridGeometry::faceIdx()][scvf.dofIndex()][0] = initSol[Indices::velocity(scvf.directionIndex())];
}
/*!
* \brief An additional drag term can be included as source term for the momentum balance
* to mimic 3D flow behavior in 2D:
* \f[
* f_{drag} = -(8 \mu / h^2)v
* \f]
* Here, \f$h\f$ corresponds to the extruded height that is
* bounded by the imaginary walls. See Flekkoy et al. (1995) \cite flekkoy1995a<BR>
* A value of 8.0 is used as a default factor, corresponding
* to the velocity profile at the center plane
* of the virtual height (maximum velocity). Setting this value to 12.0 corresponds
* to an depth-averaged velocity (Venturoli and Boek, 2006) \cite venturoli2006a.
*/
Scalar pseudo3DWallFriction(const Scalar velocity,
const Scalar viscosity,
const Scalar height,
const Scalar factor = 8.0) const
{
static_assert(dim == 2, "Pseudo 3D wall friction may only be used in 2D");
return -factor * velocity * viscosity / (height*height);
}
//! Convenience function for staggered grid implementation.
template <class ElementVolumeVariables, class ElementFaceVariables, class G = FVGridGeometry>
typename std::enable_if<G::discMethod == DiscretizationMethod::staggered, Scalar>::type
pseudo3DWallFriction(const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const Scalar height,
const Scalar factor = 8.0) const
{
const Scalar velocity = elemFaceVars[scvf].velocitySelf();
const Scalar viscosity = elemVolVars[scvf.insideScvIdx()].effectiveViscosity();
return pseudo3DWallFriction(velocity, viscosity, height, factor);
}
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
......
......@@ -145,10 +145,10 @@ public:
FacePrimaryVariables computeStorageForFace(const Problem& problem,
const SubControlVolumeFace& scvf,
const VolumeVariables& volVars,
const ElementFaceVariables& elementFaceVars) const
const ElementFaceVariables& elemFaceVars) const
{
FacePrimaryVariables storage(0.0);
const Scalar velocity = elementFaceVars[scvf].velocitySelf();
const Scalar velocity = elemFaceVars[scvf].velocitySelf();
storage[0] = volVars.density() * velocity;
return storage;
}
......@@ -159,13 +159,12 @@ public:
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elementFaceVars) const
const ElementFaceVariables& elemFaceVars) const
{
FacePrimaryVariables source(0.0);
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideVolVars = elemVolVars[insideScvIdx];
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
source += problem.gravity()[scvf.directionIndex()] * insideVolVars.density();
source += problem.source(element, fvGeometry, elemVolVars, elementFaceVars, scvf)[Indices::velocity(scvf.directionIndex())];
source += problem.source(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())];
return source;
}
......@@ -176,11 +175,11 @@ public:
const SubControlVolumeFace& scvf,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elementFaceVars,
const ElementFaceVariables& elemFaceVars,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
FluxVariables fluxVars;
return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elementFaceVars);
return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars);
}
/*!
......@@ -264,7 +263,7 @@ protected:
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elementFaceVars,
const ElementFaceVariables& elemFaceVars,
const ElementBoundaryTypes& elemBcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
......@@ -276,7 +275,7 @@ protected:
// set a fixed value for the velocity for Dirichlet boundary conditions
if(bcTypes.isDirichlet(Indices::velocity(scvf.directionIndex())))
{
const Scalar velocity = elementFaceVars[scvf].velocitySelf();
const Scalar velocity = elemFaceVars[scvf].velocitySelf();
const Scalar dirichletValue = problem.dirichlet(element, scvf)[Indices::velocity(scvf.directionIndex())];
residual = velocity - dirichletValue;
}
......@@ -286,19 +285,19 @@ protected:
{
// set the given Neumann flux for the face on the boundary itself
const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
residual = problem.neumann(element, fvGeometry, elemVolVars, elementFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
residual = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
* extrusionFactor * scvf.area();
// treat the remaining (normal) faces of the staggered control volume
FluxVariables fluxVars;
residual += fluxVars.computeNormalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elementFaceVars);
residual += fluxVars.computeNormalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars);
}
// For symmetry boundary conditions, there is no flow accross the boundary and
// we therefore treat it like a Dirichlet boundary conditions with zero velocity
if(bcTypes.isSymmetry())
{
const Scalar velocity = elementFaceVars[scvf].velocitySelf();
const Scalar velocity = elemFaceVars[scvf].velocitySelf();
const Scalar fixedValue = 0.0;
residual = velocity - fixedValue;
}
......@@ -307,7 +306,7 @@ protected:
if(bcTypes.isOutflow(Indices::velocity(scvf.directionIndex())))
{
if(bcTypes.isDirichlet(Indices::conti0EqIdx))
residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elementFaceVars, elemFluxVarsCache);
residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
else
DUNE_THROW(Dune::InvalidStateException, "Face at " << scvf.center() << " has an outflow BC for the momentum balance but no Dirichlet BC for the pressure!");
}
......
......@@ -97,3 +97,25 @@ dune_add_test(NAME test_navierstokes_angeli
--files ${CMAKE_SOURCE_DIR}/test/references/test_angeli-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_angeli-00045.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_navierstokes_angeli")
dune_add_test(NAME test_stokes_channel_3d
SOURCES test_stokes_channel_3d.cc
COMPILE_DEFINITIONS DIM_3D=1
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/stokes_channel_3d-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_stokes_channel_3d-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_stokes_channel_3d"
--zeroThreshold {"velocity_component \(m/s\)":1e-12})
dune_add_test(NAME test_stokes_channel_pseudo3d
SOURCES test_stokes_channel_3d.cc
COMPILE_DEFINITIONS DIM_3D=0
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/stokes_channel_pseudo3d-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_stokes_channel_pseudo3d-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_stokes_channel_pseudo3d"
--zeroThreshold {"velocity_component \(m/s\)":1e-12})
// -*- 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 2 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
*
* \brief Channel flow test for the staggered grid (Navier-)Stokes model.
*
* The channel is either modeled in 3D or in 2D, using an additional wall friction term
* to mimic the 3D behavior of the flow.
*
*/
#ifndef DUMUX_3D_CHANNEL_PROBLEM_HH
#define DUMUX_3D_CHANNEL_PROBLEM_HH
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/navierstokes/model.hh>
#include <dune/common/float_cmp.hh>
namespace Dumux
{
template <class TypeTag>
class ThreeDChannelTestProblem;
namespace Properties
{
NEW_TYPE_TAG(ThreeDChannelTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
// the fluid system
SET_PROP(ThreeDChannelTestTypeTag, FluidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
};
// Set the grid type
#if DIM_3D
SET_TYPE_PROP(ThreeDChannelTestTypeTag, Grid, Dune::YaspGrid<3>);
#else
SET_TYPE_PROP(ThreeDChannelTestTypeTag, Grid, Dune::YaspGrid<2>);
#endif
// Set the problem property
SET_TYPE_PROP(ThreeDChannelTestTypeTag, Problem, ThreeDChannelTestProblem<TypeTag> );
SET_BOOL_PROP(ThreeDChannelTestTypeTag, EnableFVGridGeometryCache, true);
SET_BOOL_PROP(ThreeDChannelTestTypeTag, EnableGridFluxVariablesCache, true);
SET_BOOL_PROP(ThreeDChannelTestTypeTag, EnableGridVolumeVariablesCache, true);
}
/*!
* \brief Test problem for the one-phase model:
\todo doc me!
*/
template <class TypeTag>
class ThreeDChannelTestProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using Element = typename GridView::template Codim<0>::Entity;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
static constexpr bool enablePseudoThreeDWallFriction = !DIM_3D;
public:
ThreeDChannelTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry), eps_(1e-6)
{
deltaP_ = getParam<Scalar>("Problem.DeltaP");
height_ = getParam<Scalar>("Problem.Height");
rho_ = getParam<Scalar>("Component.LiquidDensity");
nu_ = getParam<Scalar>("Component.LiquidKinematicViscosity");
if(dim == 3 && !Dune::FloatCmp::eq(height_, this->fvGridGeometry().bBoxMax()[2]))
DUNE_THROW(Dune::InvalidStateException, "z-dimension must equal height");
if(enablePseudoThreeDWallFriction)
extrusionFactor_ = 2.0/3.0 * height_;
else
extrusionFactor_ = 1.0;
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief Return the temperature within the domain in [K].
*
* This problem assumes a temperature of 10 degrees Celsius.
*/
Scalar temperature() const
{ return 273.15 + 10; } // 10C
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume face.
*/
using ParentType::source;
template<class ElementVolumeVariables, class ElementFaceVariables>
NumEqVector source(const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace &scvf) const
{
auto source = NumEqVector(0.0);
#if !DIM_3D
static const Scalar height = getParam<Scalar>("Problem.Height");
static const Scalar factor = getParam<Scalar>("Problem.PseudoWallFractionFactor", 8.0);
source[scvf.directionIndex()] = this->pseudo3DWallFriction(scvf, elemVolVars, elemFaceVars, height, factor);
#endif
return source;
}
Scalar extrusionFactorAtPos(const GlobalPosition& pos) const
{ return extrusionFactor_; }
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary control volume.
*
* \param globalPos The position of the center of the finite volume
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
// set Dirichlet values for the velocity everywhere
values.setDirichlet(Indices::momentumXBalanceIdx);
values.setDirichlet(Indices::momentumYBalanceIdx);
if(dim == 3)
values.setDirichlet(Indices::momentumZBalanceIdx);
// set a fixed pressure in one cell
if (isOutlet_(globalPos) || isInlet_(globalPos))
{
values.setDirichlet(Indices::conti0EqIdx);
values.setOutflow(Indices::momentumXBalanceIdx);
values.setOutflow(Indices::momentumYBalanceIdx);
if(dim == 3)
values.setOutflow(Indices::momentumZBalanceIdx);
}
else
values.setOutflow(Indices::conti0EqIdx);
return values;
}
/*!
* \brief Evaluate the boundary conditions for a dirichlet
* control volume.
*
* \param globalPos The center of the finite volume which ought to be set.
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
if(isInlet_(globalPos))
values[Indices::pressureIdx] = 1e5 + deltaP_;
if(isOutlet_(globalPos))
values[Indices::pressureIdx] = 1e5;
return values;
}
// \}
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param globalPos The global position
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
return values;
}
//! Returns the analytical solution for the flux through the rectangular channel
Scalar analyticalFlux() const
{
const Scalar h = height_;
const Scalar w = this->fvGridGeometry().bBoxMax()[1];
const Scalar L = this->fvGridGeometry().bBoxMax()[0];
const Scalar mu = nu_*rho_;
return h*h*h * w * deltaP_ / (12*mu*L) * (1.0 - 0.630 * h/w);
}
private:
bool isInlet_(const GlobalPosition& globalPos) const
{
return globalPos[0] < eps_;
}
bool isOutlet_(const GlobalPosition& globalPos) const
{
return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_;
}
Scalar eps_;
Scalar deltaP_;
Scalar extrusionFactor_;
Scalar height_;
Scalar rho_;
Scalar nu_;
};
} //end namespace
#endif
// -*- 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 2 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/>. *
*****************************************************************************/
#include <config.h>
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.hh>
#include "stokeschannel3dproblem.hh"
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/staggeredfvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/freeflow/navierstokes/staggered/fluxoversurface.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
*
* \param progName The name of the program, that was tried to be started.
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
*/
void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
std::string errorMessageOut = "\nUsage: ";
errorMessageOut += progName;
errorMessageOut += " [options]\n";
errorMessageOut += errorMsg;
errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
"\t-TimeManager.TEnd End of the simulation [s] \n"
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.File Name of the file containing the grid \n"
"\t definition in DGF format\n"
"\t-SpatialParams.LensLowerLeftX x-coordinate of the lower left corner of the lens [m] \n"
"\t-SpatialParams.LensLowerLeftY y-coordinate of the lower left corner of the lens [m] \n"
"\t-SpatialParams.LensUpperRightX x-coordinate of the upper right corner of the lens [m] \n"
"\t-SpatialParams.LensUpperRightY y-coordinate of the upper right corner of the lens [m] \n"
"\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
"\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
std::cout << errorMessageOut
<< "\n";
}
}
int main(int argc, char** argv) try
{
using namespace Dumux;
// define the type tag for this problem
using TypeTag = TTAG(ThreeDChannelTestTypeTag);
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
// print dumux start message
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
// parse command line arguments and input file
Parameters::init(argc, argv, usage);
// try to create a grid (from the given grid file or the input file)
using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
GridCreator::makeGrid();
GridCreator::loadBalance();
////////////////////////////////////////////////////////////
// run instationary non-linear problem on this grid
////////////////////////////////////////////////////////////
// we compute on the leaf grid view
const auto& leafGridView = GridCreator::grid().leafGridView();