// -*- 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
* \ingroup NavierStokesNCTests
* \brief Channel flow test for the multi-component staggered grid (Navier-)Stokes model
*/
#ifndef DUMUX_CHANNEL_NC_TEST_PROBLEM_HH
#define DUMUX_CHANNEL_NC_TEST_PROBLEM_HH
#ifndef ENABLECACHING
#define ENABLECACHING 1
#endif
#include <dune/grid/yaspgrid.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/h2oair.hh>
#include <dumux/material/fluidsystems/1padapter.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
namespace Dumux {
template <class TypeTag>
class ChannelNCTestProblem;
namespace Properties {
// Create new type tags
namespace TTag {
#if !NONISOTHERMAL
struct ChannelNCTest { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
#else
struct ChannelNCTest { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
#endif
} // end namespace TTag
// Select the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::ChannelNCTest>
{
using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
static constexpr int phaseIdx = H2OAir::liquidPhaseIdx;
using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
};
template<class TypeTag>
struct ReplaceCompEqIdx<TypeTag, TTag::ChannelNCTest> { static constexpr int value = 0; };
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::ChannelNCTest> { using type = Dune::YaspGrid<2>; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::ChannelNCTest> { using type = Dumux::ChannelNCTestProblem<TypeTag> ; };
template<class TypeTag>
struct EnableFVGridGeometryCache<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = ENABLECACHING; };
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = ENABLECACHING; };
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = ENABLECACHING; };
template<class TypeTag>
struct EnableGridFaceVariablesCache<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = ENABLECACHING; };
// Use mole fraction formulation
#if USE_MASS
template<class TypeTag>
struct UseMoles<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = false; };
#else
template<class TypeTag>
struct UseMoles<TypeTag, TTag::ChannelNCTest> { static constexpr bool value = true; };
#endif
}
/*!
* \ingroup NavierStokesNCTests
* \brief Test problem for the one-phase model.
* \todo doc me!
*/
template <class TypeTag>
class ChannelNCTestProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
static constexpr auto dimWorld = GetPropType<TypeTag, Properties::GridView>::dimensionworld;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
static constexpr auto compIdx = 1;
static constexpr auto transportCompIdx = Indices::conti0EqIdx + compIdx;
static constexpr auto transportEqIdx = Indices::conti0EqIdx + compIdx;
public:
ChannelNCTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry), eps_(1e-6)
{
inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
FluidSystem::init();
deltaP_.resize(this->fvGridGeometry().numCellCenterDofs());
}
/*!
* \name Problem parameters
*/
// \{
bool shouldWriteRestartFile() const
{
return false;
}
/*!
* \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 Return the sources within the domain.
*
* \param globalPos The global position
*/
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{
return NumEqVector(0.0);
}
// \}
/*!
* \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;
if(isInlet(globalPos))
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
values.setDirichlet(transportCompIdx);
#if NONISOTHERMAL
values.setDirichlet(Indices::temperatureIdx);
#endif
}
else if(isOutlet(globalPos))
{
values.setDirichlet(Indices::pressureIdx);
values.setOutflow(transportEqIdx);
#if NONISOTHERMAL
values.setOutflow(Indices::energyEqIdx);
#endif
}
else
{
// set Dirichlet values for the velocity everywhere
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
values.setNeumann(Indices::conti0EqIdx);
values.setNeumann(transportEqIdx);
#if NONISOTHERMAL
values.setNeumann(Indices::energyEqIdx);
#endif
}
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 = initialAtPos(globalPos);
// give the system some time so that the pressure can equilibrate, then start the injection of the tracer
if(isInlet(globalPos))
{
if(time() >= 10.0 || inletVelocity_ < eps_)
{
Scalar moleFracTransportedComp = 1e-3;
#if USE_MASS
Scalar averageMolarMassPhase = moleFracTransportedComp * FluidSystem::molarMass(compIdx)
+ (1. - moleFracTransportedComp) * FluidSystem::molarMass(1-compIdx);
values[transportCompIdx] = moleFracTransportedComp * FluidSystem::molarMass(compIdx)
/ averageMolarMassPhase;
#else
values[transportCompIdx] = moleFracTransportedComp;
#endif
#if NONISOTHERMAL
values[Indices::temperatureIdx] = 293.15;
#endif
}
}
return values;
}
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param globalPos The global position
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values;
values[Indices::pressureIdx] = 1.1e+5;
values[transportCompIdx] = 0.0;
#if NONISOTHERMAL
values[Indices::temperatureIdx] = 283.15;
#endif
// parabolic velocity profile
values[Indices::velocityXIdx] = inletVelocity_*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1])*(this->fvGridGeometry().bBoxMax()[1] - globalPos[1])
/ (0.25*(this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1])*(this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]));
values[Indices::velocityYIdx] = 0.0;
return values;
}
/*!
* \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
*
* \param gridVariables The grid variables
* \param sol The solution vector
*/
template<class GridVariables, class SolutionVector>
void calculateDeltaP(const GridVariables& gridVariables, const SolutionVector& sol)
{
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
auto ccDofIdx = scv.dofIndex();
auto elemVolVars = localView(gridVariables.curGridVolVars());
elemVolVars.bindElement(element, fvGeometry, sol);
deltaP_[ccDofIdx] = elemVolVars[scv].pressure() - 1.1e5;
}
}
}
auto& getDeltaP() const
{ return deltaP_; }
// \}
void setTimeLoop(TimeLoopPtr timeLoop)
{
timeLoop_ = timeLoop;
}
Scalar time() const
{
return timeLoop_->time();
}
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_;
}
const Scalar eps_;
Scalar inletVelocity_;
TimeLoopPtr timeLoop_;
std::vector<Scalar> deltaP_;
};
} //end namespace
#endif