diff --git a/dumux/freeflow/staggeredni/fluxvariables.hh b/dumux/freeflow/staggeredni/fluxvariables.hh index ac33aceab492309b4a0c4a293ab7f9a58e81274b..661f02221ec3df51986a4218fe451c20f79dc4c3 100644 --- a/dumux/freeflow/staggeredni/fluxvariables.hh +++ b/dumux/freeflow/staggeredni/fluxvariables.hh @@ -33,24 +33,16 @@ namespace Dumux namespace Properties { // forward declaration -//NEW_PROP_TAG(EnableComponentTransport); NEW_PROP_TAG(EnableEnergyBalanceStokes); NEW_PROP_TAG(EnableInertiaTerms); } -// // forward declaration -// template<class TypeTag, bool enableComponentTransport, bool enableEnergyBalance> -// class FreeFlowFluxVariablesImpl; - /*! * \ingroup ImplicitModel * \brief The flux variables class * specializations are provided for combinations of physical processes * \note Not all specializations are currently implemented */ -// template<class TypeTag> -// using FreeFlowFluxVariables = FreeFlowFluxVariablesImpl<TypeTag, GET_PROP_VALUE(TypeTag, EnableComponentTransport), -// GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>; // specialization for immiscible, non-isothermal flow template<class TypeTag> @@ -69,7 +61,6 @@ class FreeFlowFluxVariablesImpl<TypeTag, false, true> : public FreeFlowFluxVaria using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables); using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType); - static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms); static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles); using ParentType = FreeFlowFluxVariablesImpl<TypeTag, false, false>; @@ -79,11 +70,20 @@ class FreeFlowFluxVariablesImpl<TypeTag, false, true> : public FreeFlowFluxVaria dim = GridView::dimension, dimWorld = GridView::dimensionworld, - momentumBalanceIdx = Indices::momentumBalanceIdx, - energyBalanceIdx = Indices::energyBalanceIdx, + energyBalanceIdx = Indices::energyBalanceIdx }; public: + /*! + * \brief Returns the fluxes for the cell center primary variables + * \param problem The problem + * \param element The element + * \param fvGeometry The finite-volume geometry + * \param elemVolVars All volume variables for the element + * \param globalFaceVars The face variables + * \param scvf The sub control volume face + * \param fluxVarsCache The flux variables cache + */ CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem, const Element &element, const FVElementGeometry& fvGeometry, @@ -93,13 +93,21 @@ public: const FluxVariablesCache& fluxVarsCache) { CellCenterPrimaryVariables flux(0.0); + flux = ParentType::computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, globalFaceVars, scvf, fluxVarsCache); flux += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, globalFaceVars, scvf); flux += HeatConductionType::diffusiveFluxForCellCenter(problem, element, fvGeometry, elemVolVars, scvf); return flux; } private: - + /*! + * \brief Returns the advective fluxes for the cell center primary variables + * \param problem The problem + * \param fvGeometry The finite-volume geometry + * \param elemVolVars All volume variables for the element + * \param globalFaceVars The face variables + * \param scvf The sub control volume face + */ CellCenterPrimaryVariables advectiveFluxForCellCenter_(const Problem& problem, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, @@ -117,7 +125,7 @@ private: if(scvf.boundary()) { const auto bcTypes = problem.boundaryTypesAtPos(scvf.center()); - if(bcTypes.isOutflow(momentumBalanceIdx)) + if(bcTypes.isOutflow(energyBalanceIdx)) isOutflow = true; } @@ -127,7 +135,7 @@ private: const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity); const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars; - const auto& downstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars; + const auto& downstreamVolVars = insideIsUpstream ? outsideVolVars : insideVolVars; const Scalar upWindWeight = GET_PROP_VALUE(TypeTag, ImplicitUpwindWeight); const Scalar upstreamDensity = useMoles ? upstreamVolVars.molarDensity() : upstreamVolVars.density(); @@ -138,8 +146,8 @@ private: flux[energyBalanceIdx] = (upWindWeight * upstreamDensity * upstreamEnthalpy + (1.0 - upWindWeight) * downstreamDensity * downstreamEnthalpy) * velocity; - flux *= scvf.area() * sign(scvf.outerNormalScalar()); + return flux; } }; diff --git a/dumux/freeflow/staggeredni/localresidual.hh b/dumux/freeflow/staggeredni/localresidual.hh index fb3276b461bcbede296898dc6ee2d3167be979da..43b22fe26249ab5bf89baed2330bcc411271af7c 100644 --- a/dumux/freeflow/staggeredni/localresidual.hh +++ b/dumux/freeflow/staggeredni/localresidual.hh @@ -20,7 +20,7 @@ * \file * * \brief Element-wise calculation of the Jacobian matrix for problems - * using the non-isothermal Stokes box model with a staggered grid. + * using the non-isothermal Stokes model with a staggered grid. * */ #ifndef DUMUX_STAGGERED_NAVIERSTOKES_NI_LOCAL_RESIDUAL_HH @@ -39,10 +39,8 @@ namespace Dumux namespace Properties { // forward declaration -NEW_PROP_TAG(EnableComponentTransport); NEW_PROP_TAG(EnableEnergyBalanceStokes); NEW_PROP_TAG(EnableInertiaTerms); -NEW_PROP_TAG(ReplaceCompEqIdx); } /*! @@ -83,29 +81,63 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false, true> : public Staggered using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars); - static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Element = typename GridView::template Codim<0>::Entity; + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); + +// static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles); public: /*! - * \brief Evaluate the amount the additional quantities to the Stokes model - * (energy equation). + * \brief Evaluate the cell center storage terms of the non-isothermal Stokes model. * * The result should be averaged over the volume (e.g. phase mass * inside a sub control volume divided by the volume) + * + * \param scv The sub control volume + * \param volVars The volume variables */ CellCenterPrimaryVariables computeStorageForCellCenter(const SubControlVolume& scv, const VolumeVariables& volVars) { CellCenterPrimaryVariables storage(0.0); - const Scalar density = useMoles? volVars.molarDensity() : volVars.density(); +// const Scalar density = useMoles? volVars.molarDensity() : volVars.density(); // compute storage of mass - storage[massBalanceIdx] = volVars.density(); + storage = ParentType::computeStorageForCellCenter(scv, volVars); // compute the storage of energy - storage[energyBalanceIdx] = density * volVars.internalEnergy(); + storage[energyBalanceIdx] = volVars.density() * volVars.internalEnergy(); return storage; } + + + /*! + * \brief Evaluate the cell center source terms of the non-isothermal Stokes model. + * + * The result should be averaged over the volume (e.g. phase mass + * inside a sub control volume divided by the volume) + * + * \param element The element + * \param fvGeometry The finite-volume geometry + * \param elemVolVars All volume variables for the element + * \param globalFaceVars The face variables + * \param scv The sub control volume + */ + CellCenterPrimaryVariables computeSourceForCellCenter(const Element &element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const GlobalFaceVars& globalFaceVars, + const SubControlVolume &scv) + { + CellCenterPrimaryVariables source(0.0); + + source = this->problem().sourceAtPos(scv.center())[cellCenterIdx][massBalanceIdx]; + source = this->problem().sourceAtPos(scv.center())[cellCenterIdx][energyBalanceIdx]; + + return source; + } }; } // end namespace diff --git a/dumux/freeflow/staggeredni/model.hh b/dumux/freeflow/staggeredni/model.hh index 0f8b3a7a3a41751779313e4205cce206d7502df7..583226c8ca46495e535290c1bb58d5d5b83fa44d 100644 --- a/dumux/freeflow/staggeredni/model.hh +++ b/dumux/freeflow/staggeredni/model.hh @@ -37,20 +37,7 @@ namespace Dumux * \ingroup NavierStokesModel * \brief A single-phase, non-isothermal flow model using the fully implicit scheme. * - * Single-phase, non-isothermal flow model, which uses a standard Darcy approach as the - * equation for the conservation of momentum: - * \f[ - v = - \frac{\textbf K}{\mu} - \left(\textbf{grad}\, p - \varrho {\textbf g} \right) - * \f] - * - * and solves the mass continuity equation: - * \f[ - \phi \frac{\partial \varrho}{\partial t} + \text{div} \left\lbrace - - \varrho \frac{\textbf K}{\mu} \left( \textbf{grad}\, p -\varrho {\textbf g} \right) \right\rbrace = q, - * \f] - * All equations are discretized using a vertex-centered finite volume (box) - * or cell-centered finite volume scheme as spatial + * All equations are discretized using a staggered grid as spatial * and the implicit Euler method as time discretization. * The model supports compressible as well as incompressible fluids. */ @@ -95,7 +82,6 @@ public: // add temperature to output auto& vtkOutputModule = problem.vtkOutputModule(); -// vtkOutputModule.addPrimaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature();}); vtkOutputModule.addPrimaryVariable("temperature", Indices::temperatureIdx); } }; diff --git a/dumux/freeflow/staggeredni/properties.hh b/dumux/freeflow/staggeredni/properties.hh index dbe0b517b5d2dc8c804066b3c0b64ffa705323aa..19f6a6bf2ef455018952cb967440e074c999f239 100644 --- a/dumux/freeflow/staggeredni/properties.hh +++ b/dumux/freeflow/staggeredni/properties.hh @@ -28,7 +28,6 @@ #define DUMUX_NAVIERSTOKES_NI_PROPERTIES_HH #include <dumux/freeflow/staggered/properties.hh> -//#include <dumux/porousmediumflow/nonisothermal/implicit/properties.hh> namespace Dumux { @@ -53,22 +52,6 @@ NEW_PROP_TAG(EnableEnergyBalanceStokes); ////////////////////////////////////////////////////////////////// // Property tags ////////////////////////////////////////////////////////////////// - -//NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system -//NEW_PROP_TAG(Indices); //!< Enumerations for the model -//NEW_PROP_TAG(FluidSystem); //!< The type of the fluid system to use -//NEW_PROP_TAG(Fluid); //!< The fluid used for the default fluid system -//NEW_PROP_TAG(FluidState); //!< The type of the fluid state to use -//NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem -//NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< Returns weight of the upwind cell when calculating fluxes -//NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation -//NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output -//NEW_PROP_TAG(EnableInertiaTerms); //!< Returns whether to include inertia terms in the momentum balance eq or not (Stokes / Navier-Stokes) -//NEW_PROP_TAG(BoundaryValues); //!< Type to set values on the boundary -//NEW_PROP_TAG(EnableComponentTransport); //!< Returns whether to consider component transport or not -//NEW_PROP_TAG(EnableEnergyTransport); //!< Returns whether to consider energy transport or not -//NEW_PROP_TAG(FaceVariables); //!< Returns whether to consider energy transport or not -//NEW_PROP_TAG(UseMoles); //!< Defines whether molar (true) or mass (false) density is used NEW_PROP_TAG(PhaseIdx); //!< Defines the phaseIdx // \} diff --git a/dumux/freeflow/staggeredni/propertydefaults.hh b/dumux/freeflow/staggeredni/propertydefaults.hh index dc1917b603b7b4c596235217dfb84c05c5321e0b..93163c2866e2c24754072d742dec04694220d561 100644 --- a/dumux/freeflow/staggeredni/propertydefaults.hh +++ b/dumux/freeflow/staggeredni/propertydefaults.hh @@ -39,10 +39,10 @@ #include "../staggered/propertydefaults.hh" #include <dumux/implicit/staggered/localresidual.hh> -#include <dumux/material/fluidsystems/gasphase.hh> +//#include <dumux/material/fluidsystems/gasphase.hh> #include <dumux/material/fluidsystems/liquidphase.hh> -#include <dumux/material/components/nullcomponent.hh> -#include <dumux/material/fluidsystems/1p.hh> +//#include <dumux/material/components/nullcomponent.hh> +//#include <dumux/material/fluidsystems/1p.hh> #include <dumux/material/fluidstates/immiscible.hh> @@ -78,110 +78,6 @@ SET_TYPE_PROP(NavierStokesNI, HeatConductionType, FouriersLaw<TypeTag>); SET_INT_PROP(NavierStokesNI, PhaseIdx, 0); //!< Defines the phaseIdx - -////! average is used as default model to compute the effective thermal heat conductivity -//SET_PROP(NavierStokesNI, ThermalConductivityModel) -//{ private : -// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; -// public: -// typedef ThermalConductivityAverage<Scalar> type; -//}; - - -/*! - * \brief The fluid state which is used by the volume variables to - * store the thermodynamic state. This should be chosen - * appropriately for the model ((non-)isothermal, equilibrium, ...). - * This can be done in the problem. - */ -//SET_PROP(NavierStokesNI, FluidState) -//{ -// private: -// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; -// typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; -// public: -// typedef ImmiscibleFluidState<Scalar, FluidSystem> type; -//}; - -// //! Enable advection -// SET_BOOL_PROP(NavierStokes, EnableAdvection, true); -// -// //! The one-phase model has no molecular diffusion -// SET_BOOL_PROP(NavierStokes, EnableMolecularDiffusion, false); -// -//! Non-Isothermal model by default -//SET_BOOL_PROP(NavierStokesNI, EnableEnergyBalance, true); -// -// //! The indices required by the isothermal single-phase model -// SET_TYPE_PROP(NavierStokes, Indices, NavierStokesCommonIndices<TypeTag>); -// -// //! The weight of the upwind control volume when calculating -// //! fluxes. Use central differences by default. -// SET_SCALAR_PROP(NavierStokes, ImplicitMassUpwindWeight, 0.5); -// -// //! weight for the upwind mobility in the velocity calculation -// //! fluxes. Use central differences by default. -// SET_SCALAR_PROP(NavierStokes, ImplicitMobilityUpwindWeight, 0.5); -// -// //! The fluid system to use by default -// SET_TYPE_PROP(NavierStokes, FluidSystem, Dumux::FluidSystems::OneP<typename GET_PROP_TYPE(TypeTag, Scalar), typename GET_PROP_TYPE(TypeTag, Fluid)>); -// -// SET_PROP(NavierStokes, Fluid) -// { private: -// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; -// public: -// typedef FluidSystems::LiquidPhase<Scalar, Dumux::NullComponent<Scalar> > type; -// }; -// -// /*! -// * \brief The fluid state which is used by the volume variables to -// * store the thermodynamic state. This should be chosen -// * appropriately for the model ((non-)isothermal, equilibrium, ...). -// * This can be done in the problem. -// */ -// SET_PROP(NavierStokes, FluidState){ -// private: -// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; -// typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; -// public: -// typedef Dumux::ImmiscibleFluidState<Scalar, FluidSystem> type; -// }; -// -// // disable velocity output by default -// SET_BOOL_PROP(NavierStokes, VtkAddVelocity, true); -// -// // enable gravity by default -// SET_BOOL_PROP(NavierStokes, ProblemEnableGravity, true); -// -// SET_BOOL_PROP(NavierStokes, EnableInertiaTerms, true); -// -// SET_BOOL_PROP(NavierStokes, EnableEnergyTransport, true); -// - -////////////////////////////////////////////////////////////////// -// Property values for isothermal model required for the general non-isothermal model -////////////////////////////////////////////////////////////////// - -// set isothermal Model -// SET_TYPE_PROP(NavierStokesNI, IsothermalModel, NavierStokesModel<TypeTag>); - -//set isothermal VolumeVariables -// SET_TYPE_PROP(NavierStokesNI, IsothermalVolumeVariables, NavierStokesVolumeVariables<TypeTag>); - -//set isothermal LocalResidual -// SET_TYPE_PROP(NavierStokesNI, IsothermalLocalResidual, ImmiscibleLocalResidual<TypeTag>); - -//set isothermal Indices -// SET_TYPE_PROP(NavierStokesNI, IsothermalIndices, NavierStokesCommonIndices<TypeTag>); - -//set isothermal NumEq -// SET_INT_PROP(NavierStokesNI, IsothermalNumEq, 1); - -//set non-isothermal NumEq -// SET_INT_PROP(NavierStokesNI, NonIsothermalNumEq, 1); - - -// \} } // end namespace Properties } // end namespace Dumux diff --git a/dumux/freeflow/staggeredni/volumevariables.hh b/dumux/freeflow/staggeredni/volumevariables.hh index d87b0d272284461f9845b80dd44106b5007760c2..3abb65433e5e25893fffe322ca5545cd2a3aeeee 100644 --- a/dumux/freeflow/staggeredni/volumevariables.hh +++ b/dumux/freeflow/staggeredni/volumevariables.hh @@ -73,7 +73,7 @@ public: { ParentType::update(elemSol, problem, element, scv); - completeFluidState(elemSol, problem, element, scv, fluidState_); + completeFluidState(elemSol, problem, element, scv, this->fluidState_); } /*! @@ -86,25 +86,25 @@ public: FluidState& fluidState) { fluidState.setTemperature(elemSol[0][Indices::temperatureIdx]); - fluidState.setPressure(/*phaseIdx=*/0, elemSol[0][Indices::pressureIdx]); + fluidState.setPressure(phaseIdx, elemSol[0][Indices::pressureIdx]); // saturation in a single phase is always 1 and thus redundant // to set. But since we use the fluid state shared by the // immiscible multi-phase models, so we have to set it here... - fluidState.setSaturation(/*phaseIdx=*/0, 1.0); + fluidState.setSaturation(phaseIdx, 1.0); typename FluidSystem::ParameterCache paramCache; - paramCache.updatePhase(fluidState, /*phaseIdx=*/0); + paramCache.updatePhase(fluidState, phaseIdx); - Scalar value = FluidSystem::density(fluidState, paramCache, /*phaseIdx=*/0); - fluidState.setDensity(/*phaseIdx=*/0, value); + Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx); + fluidState.setDensity(phaseIdx, value); - value = FluidSystem::viscosity(fluidState, paramCache, /*phaseIdx=*/0); - fluidState.setViscosity(/*phaseIdx=*/0, value); + value = FluidSystem::viscosity(fluidState, paramCache, phaseIdx); + fluidState.setViscosity(phaseIdx, value); // compute and set the enthalpy - value = FluidSystem::enthalpy(fluidState, paramCache, /*phaseIdx=*/0); - fluidState.setEnthalpy(/*phaseIdx=*/0, value); + value = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx); + fluidState.setEnthalpy(phaseIdx, value); } /*! @@ -125,6 +125,17 @@ public: return this->fluidState_.molarDensity(phaseIdx); } + /*! + * \brief Returns the molar mass of a given phase within the + * control volume. + * + * \param phaseIdx The phase index + */ + Scalar molarMass() const + { + return this->fluidState_.averageMolarMass(phaseIdx); + } + /*! * \brief Returns the total internal energy of the fluid phase in the * sub-control volume. @@ -153,30 +164,19 @@ public: Scalar thermalConductivity() const { return FluidSystem::thermalConductivity(this->fluidState_, phaseIdx); } + /*! + * \brief Returns the temperature \f$\mathrm{[K]}\f$ + * of the fluid phase in the sub-control volume. + */ + Scalar temperature() const + { return this->fluidState_.temperature(phaseIdx); } protected: - FluidState fluidState_; // this method gets called by the parent class. since this method // is protected, we are friends with our parent... friend class NavierStokesVolumeVariables<TypeTag>; - static Scalar temperature_(const CellCenterPrimaryVariables &priVars, - const Problem& problem, - const Element &element, - const FVElementGeometry &fvGeometry, - const int scvIdx) - { - return priVars[temperatureIdx]; - } - - template<class ParameterCache> - static Scalar enthalpy_(const FluidState& fluidState, - const ParameterCache& paramCache, - const int phaseIdx) - { - return FluidSystem::enthalpy(fluidState, paramCache, phaseIdx); - } }; } // end namespace