Commit f3ecfd37 authored by Alexander Kissinger's avatar Alexander Kissinger
Browse files

- Added the general non-isothermal model in implicit/nonisothermal

- The implicit 1p2c model is ready to use with the non-isothermal model
- Added 1p fluidmatrixinteractions to calculate effective thermal conductivity of fluid and solid
- Adjusted test/implicit/1p2c:
	- 1p2coutflowproblem.hh can be switched to nonisothermal
	- 1p2cconductionproblem.hh compares the model results to a simple analytical model
	- 1p2cconvectionproblem.hh compares the model results to another simple analytical model

Reviewed by Bernd

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@13264 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent cb7521f3
......@@ -60,6 +60,7 @@ AC_CONFIG_FILES([dumux.pc
dumux/implicit/common/Makefile
dumux/implicit/co2/Makefile
dumux/implicit/co2ni/Makefile
dumux/implicit/nonisothermal/Makefile
dumux/implicit/mpnc/Makefile
dumux/implicit/mpnc/diffusion/Makefile
dumux/implicit/mpnc/energy/Makefile
......@@ -71,7 +72,8 @@ AC_CONFIG_FILES([dumux.pc
dumux/material/binarycoefficients/Makefile
dumux/material/components/Makefile
dumux/material/components/iapws/Makefile
dumux/material/fluidmatrixinteractions/Makefile
dumux/material/fluidmatrixinteractions/Makefile
dumux/material/fluidmatrixinteractions/1p/Makefile
dumux/material/fluidmatrixinteractions/2p/Makefile
dumux/material/fluidmatrixinteractions/2pia/Makefile
dumux/material/fluidmatrixinteractions/3p/Makefile
......
......@@ -246,6 +246,20 @@ public:
int downstreamIdx() const
{ return downstreamIdx_; }
/*!
* \brief Return the local index of the upstream control volume
* for a given phase.
*/
int upstreamIdx(int phaseIdx) const
{ return upstreamIdx_; }
/*!
* \brief Return the local index of the downstream control volume
* for a given phase.
*/
int downstreamIdx(int phaseIdx) const
{ return downstreamIdx_; }
/*!
* \brief Return the volumetric flux over a face of a given phase.
*
......
......@@ -47,6 +47,11 @@ namespace Properties
NEW_TYPE_TAG(OnePTwoC);
NEW_TYPE_TAG(BoxOnePTwoC, INHERITS_FROM(BoxModel, OnePTwoC));
NEW_TYPE_TAG(CCOnePTwoC, INHERITS_FROM(CCModel, OnePTwoC));
//! The type tags for the implicit one-phase two-component non-isothermal problems
NEW_TYPE_TAG(NonIsothermal, INHERITS_FROM(OnePTwoC));
NEW_TYPE_TAG(OnePTwoCNI, INHERITS_FROM(NonIsothermal));
NEW_TYPE_TAG(BoxOnePTwoCNI, INHERITS_FROM(BoxModel, OnePTwoCNI));
NEW_TYPE_TAG(CCOnePTwoCNI, INHERITS_FROM(CCModel, OnePTwoCNI));
//////////////////////////////////////////////////////////////////
// Property tags
......@@ -66,6 +71,19 @@ NEW_PROP_TAG(UseMoles); //!< Defines whether mole (true) or mass (false) fractio
NEW_PROP_TAG(Scaling); //!< Defines Scaling of the model
NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
//////////////////////////////////////////////////////////////////
// Property tags required for the non-isothermal model
//////////////////////////////////////////////////////////////////
NEW_PROP_TAG(IsothermalModel);
NEW_PROP_TAG(IsothermalFluxVariables);
NEW_PROP_TAG(IsothermalVolumeVariables);
NEW_PROP_TAG(IsothermalLocalResidual);
NEW_PROP_TAG(IsothermalIndices);
NEW_PROP_TAG(IsothermalNumEq);
NEW_PROP_TAG(HaveVariableFormulation);
NEW_PROP_TAG(ThermalConductivityModel); //!< The model for the effective thermal conductivity
}
// \}
}
......
......@@ -37,8 +37,10 @@
#include "1p2cfluxvariables.hh"
#include "1p2cindices.hh"
#include <dumux/implicit/nonisothermal/nipropertydefaults.hh>
#include <dumux/material/spatialparams/implicitspatialparams1p.hh>
#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
#include <dumux/material/fluidmatrixinteractions/1p/thermalconductivityaverage.hh>
namespace Dumux
{
......@@ -49,23 +51,24 @@ namespace Properties
// Property values
//////////////////////////////////////////////////////////////////
SET_INT_PROP(OnePTwoC, NumEq, 2); //!< set the number of equations to 2
SET_INT_PROP(OnePTwoC, NumEq, GET_PROP_VALUE(TypeTag, IsothermalNumEq)); //!< set the number of equations to 2
SET_INT_PROP(OnePTwoC, NumPhases, 1); //!< The number of phases in the 1p2c model is 1
SET_INT_PROP(OnePTwoC, NumComponents, 2); //!< The number of components in the 1p2c model is 2
SET_SCALAR_PROP(OnePTwoC, Scaling, 1); //!< Scaling of the model is set to 1 by default
SET_BOOL_PROP(OnePTwoC, UseMoles, true); //!< Define that mole fractions are used in the balance equations
//! Use the 1p2c local residual function for the 1p2c model
SET_TYPE_PROP(OnePTwoC, LocalResidual, OnePTwoCLocalResidual<TypeTag>);
SET_TYPE_PROP(OnePTwoC, LocalResidual, typename GET_PROP_TYPE(TypeTag, IsothermalLocalResidual));
//! define the model
SET_TYPE_PROP(OnePTwoC, Model, OnePTwoCModel<TypeTag>);
SET_TYPE_PROP(OnePTwoC, Model, typename GET_PROP_TYPE(TypeTag, IsothermalModel));
//! define the VolumeVariables
SET_TYPE_PROP(OnePTwoC, VolumeVariables, OnePTwoCVolumeVariables<TypeTag>);
SET_TYPE_PROP(OnePTwoC, VolumeVariables, typename GET_PROP_TYPE(TypeTag, IsothermalVolumeVariables));
//! define the FluxVariables
SET_TYPE_PROP(OnePTwoC, FluxVariables, OnePTwoCFluxVariables<TypeTag>);
SET_TYPE_PROP(OnePTwoC, FluxVariables, typename GET_PROP_TYPE(TypeTag, IsothermalFluxVariables));
//! set default upwind weight to 1.0, i.e. fully upwind
SET_SCALAR_PROP(OnePTwoC, ImplicitMassUpwindWeight, 1.0);
......@@ -74,8 +77,7 @@ SET_SCALAR_PROP(OnePTwoC, ImplicitMassUpwindWeight, 1.0);
SET_SCALAR_PROP(OnePTwoC, ImplicitMobilityUpwindWeight, 1.0);
//! Set the indices used by the 1p2c model
SET_TYPE_PROP(OnePTwoC, Indices, Dumux::OnePTwoCIndices<TypeTag, 0>);
SET_TYPE_PROP(OnePTwoC, Indices, typename GET_PROP_TYPE(TypeTag, IsothermalIndices));
//! The spatial parameters to be employed.
//! Use ImplicitSpatialParamsOneP by default.
SET_TYPE_PROP(OnePTwoC, SpatialParams, ImplicitSpatialParamsOneP<TypeTag>);
......@@ -103,6 +105,39 @@ SET_BOOL_PROP(OnePTwoC, ProblemEnableGravity, true);
// porous medium. Taking it as a constant is only a first approximation
// (Nield, Bejan, Convection in porous media, 2006, p. 10)
SET_SCALAR_PROP(BoxModel, SpatialParamsForchCoeff, 0.55);
//! average is used as default model to compute the effective thermal heat conductivity
SET_PROP(NonIsothermal, ThermalConductivityModel)
{ private :
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
public:
typedef ThermalConductivityAverage<Scalar, VolumeVariables> type;
};
//////////////////////////////////////////////////////////////////
// Property values for isothermal model required for the general non-isothermal model
//////////////////////////////////////////////////////////////////
// set isothermal Model
SET_TYPE_PROP(OnePTwoC, IsothermalModel, OnePTwoCModel<TypeTag>);
// set isothermal FluxVariables
SET_TYPE_PROP(OnePTwoC, IsothermalFluxVariables, OnePTwoCFluxVariables<TypeTag>);
//set isothermal VolumeVariables
SET_TYPE_PROP(OnePTwoC, IsothermalVolumeVariables, OnePTwoCVolumeVariables<TypeTag>);
//set isothermal LocalResidual
SET_TYPE_PROP(OnePTwoC, IsothermalLocalResidual, OnePTwoCLocalResidual<TypeTag>);
//set isothermal Indices
SET_TYPE_PROP(OnePTwoC, IsothermalIndices, OnePTwoCIndices<TypeTag>);
//set isothermal NumEq
SET_INT_PROP(OnePTwoC, IsothermalNumEq, 2);
}
// \}
}
......
......@@ -129,9 +129,10 @@ public:
const int scvIdx,
FluidState& fluidState)
{
Scalar T = Implementation::temperature_(priVars, problem, element,
Scalar t = Implementation::temperature_(priVars, problem, element,
fvGeometry, scvIdx);
fluidState.setTemperature(T);
fluidState.setTemperature(t);
fluidState.setSaturation(phaseIdx, 1.);
fluidState.setPressure(phaseIdx, priVars[pressureIdx]);
......@@ -157,6 +158,10 @@ public:
fluidState.setDensity(phaseIdx, value);
value = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
fluidState.setViscosity(phaseIdx, value);
// compute and set the enthalpy
Scalar h = Implementation::enthalpy_(fluidState, paramCache, phaseIdx);
fluidState.setEnthalpy(phaseIdx, h);
}
/*!
......@@ -254,6 +259,14 @@ protected:
return problem.temperatureAtPos(fvGeometry.subContVol[scvIdx].global);
}
template<class ParameterCache>
static Scalar enthalpy_(const FluidState& fluidState,
const ParameterCache& paramCache,
const int phaseIdx)
{
return 0;
}
/*!
* \brief Called by update() to compute the energy related quantities.
*/
......
......@@ -14,6 +14,7 @@ SUBDIRS = 1p \
co2ni \
common \
mpnc \
nonisothermal \
richards
implicitdir = $(includedir)/dumux/implicit
......
nidir = $(includedir)/dumux/implicit/nonisothermal
ni_HEADERS = *.hh
include $(top_srcdir)/am/global-rules
// -*- 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 This file contains data which is required to calculate
* the heat fluxes over a face of a finite volume.
*
* This means temperature gradients and the normal matrix
* heat flux.
*/
#ifndef DUMUX_NI_FLUX_VARIABLES_HH
#define DUMUX_NI_FLUX_VARIABLES_HH
#include <dumux/common/math.hh>
namespace Dumux
{
/*!
* \ingroup TwoPTwoCNIModel
* \ingroup ImplicitFluxVariables
* \brief This template class contains data which is required to
* calculate the heat fluxes over a face of a finite
* volume for the non-isothermal two-phase two-component model.
* The mass fluxes are computed in the parent class.
*
* This means temperature gradients and the normal matrix
* heat flux.
*/
template <class TypeTag>
class NIFluxVariables : public GET_PROP_TYPE(TypeTag, IsothermalFluxVariables)
{
typedef typename GET_PROP_TYPE(TypeTag, IsothermalFluxVariables) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel) ThermalConductivityModel;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
enum { dimWorld = GridView::dimensionworld };
enum { dim = GridView::dimension };
typedef Dune::FieldVector<Scalar, dim> DimVector;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
// static const bool isNonIsothermal = GET_PROP_VALUE(TypeTag, IsNonIsothermal);
public:
/*!
* \brief The constructor
*
* \param problem The problem
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the fully implicit scheme
* \param faceIdx The local index of the sub-control-volume face
* \param elemVolVars The volume variables of the current element
* \param onBoundary Distinguishes if we are on a sub-control-volume face or on a boundary face
*/
NIFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const int faceIdx,
const ElementVolumeVariables &elemVolVars,
bool onBoundary = false)
: ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
{
faceIdx_ = faceIdx;
calculateValues_(problem, element, elemVolVars);
}
/*!
* \brief The total heat flux \f$\mathrm{[J/s]}\f$ due to heat conduction
* of the rock matrix over the sub-control volume face in
* direction of the face normal.
*/
Scalar normalMatrixHeatFlux() const
{ return normalMatrixHeatFlux_; }
/*!
* \brief The local temperature gradient at the IP of the considered scv face.
*/
DimVector temperatureGradient() const
{ return temperatureGrad_; }
/*!
* \brief The harmonically averaged effective thermal conductivity.
*/
Scalar effThermalConductivity() const
{ return lambdaEff_; }
protected:
void calculateValues_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
// calculate temperature gradient using finite element
// gradients
temperatureGrad_ = 0;
DimVector tmp(0.0);
for (unsigned int idx = 0; idx < this->face().numFap; idx++)
{
tmp = this->face().grad[idx];
// index for the element volume variables
int volVarsIdx = this->face().fapIndices[idx];
tmp *= elemVolVars[volVarsIdx].temperature();
temperatureGrad_ += tmp;
}
lambdaEff_ = 0;
calculateEffThermalConductivity_(problem, element, elemVolVars);
// project the heat flux vector on the face's normal vector
normalMatrixHeatFlux_ = temperatureGrad_* this->face().normal;
normalMatrixHeatFlux_ *= -lambdaEff_;
}
void calculateEffThermalConductivity_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
const unsigned i = this->face().i;
const unsigned j = this->face().j;
Scalar lambdaI, lambdaJ;
if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
{
lambdaI =
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i],
problem.spatialParams().thermalConductivitySolid(element, this->fvGeometry_, i),
problem.spatialParams().porosity(element, this->fvGeometry_, i));
lambdaJ =
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j],
problem.spatialParams().thermalConductivitySolid(element, this->fvGeometry_, j),
problem.spatialParams().porosity(element, this->fvGeometry_, j));
}
else
{
const Element& elementI = *this->fvGeometry_.neighbors[i];
FVElementGeometry fvGeometryI;
fvGeometryI.subContVol[0].global = elementI.geometry().center();
lambdaI =
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i],
problem.spatialParams().thermalConductivitySolid(elementI, fvGeometryI, 0),
problem.spatialParams().porosity(elementI, fvGeometryI, 0));
const Element& elementJ = *this->fvGeometry_.neighbors[j];
FVElementGeometry fvGeometryJ;
fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
lambdaJ =
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j],
problem.spatialParams().thermalConductivitySolid(elementJ, fvGeometryJ, 0),
problem.spatialParams().porosity(elementJ, fvGeometryJ, 0));
}
// -> harmonic mean
lambdaEff_ = harmonicMean(lambdaI, lambdaJ);
}
private:
Scalar lambdaEff_;
Scalar normalMatrixHeatFlux_;
DimVector temperatureGrad_;
int faceIdx_;
};
} // 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/>. *
*****************************************************************************/
/*!
* \file
* \brief Defines the indices used by the non-isothermal two-phase two-component model
*/
#ifndef DUMUX_NI_INDICES_HH
#define DUMUX_NI_INDICES_HH
namespace Dumux
{
/*!
* \ingroup TwoPTwoCNIModel
* \ingroup ImplicitIndices
* \brief Indices for the non-isothermal two-phase two-component model
*
* \tparam formulation The formulation, either pwsn or pnsw.
* \tparam PVOffset The first index in a primary variable vector.
*/
template <class TypeTag, int PVOffset = 0>
class NIIndices : public GET_PROP_TYPE(TypeTag, IsothermalIndices)
{
public:
static const int numEq = GET_PROP_VALUE(TypeTag, NumEq);
static const int temperatureIdx = PVOffset + numEq -1; //! The index for temperature in primary variable vectors.
static const int energyEqIdx = PVOffset + numEq -1; //! The index for energy in equation vectors.
};
}
#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/>. *
*****************************************************************************/
/*!
* \file
*
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the non-isothermal two-phase two-component fully implicit model.
*
*/
#ifndef DUMUX_NEW_NI_LOCAL_RESIDUAL_HH
#define DUMUX_NEW_NI_LOCAL_RESIDUAL_HH
namespace Dumux
{
/*!
* \ingroup TwoPTwoCNIModel
* \ingroup ImplicitLocalResidual
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the two-phase two-component fully implicit model.
*/
template<class TypeTag>
class NILocalResidual : public GET_PROP_TYPE(TypeTag, IsothermalLocalResidual)
{
typedef typename GET_PROP_TYPE(TypeTag, IsothermalLocalResidual) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
energyEqIdx = Indices::energyEqIdx,
temperatureIdx = Indices::temperatureIdx,
};
public:
/*!
* \brief Constructor. Sets the upwind weight.
*/
NILocalResidual()
{
// retrieve the upwind weight for the mass conservation equations. Use the value
// specified via the property system as default, and overwrite
// it by the run-time parameter from the Dune::ParameterTree
massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
};
/*!
* \brief Evaluate the amount of all conservation quantities
* (e.g. phase mass) within a sub-control volume.
*
* The result should be averaged over the volume (e.g. phase mass
* inside a sub control volume divided by the volume)
*
* \param storage The storage of the conservation quantity (mass or energy) within the sub-control volume
* \param scvIdx The sub-control-volume index
* \param usePrevSol Evaluate function with solution of current or previous time step
*/
void computeStorage(PrimaryVariables &storage, const int scvIdx, bool usePrevSol) const
{
// compute the storage term for phase mass
ParentType::computeStorage(storage, scvIdx, usePrevSol);
// if flag usePrevSol is set, the solution from the previous
// time step is used, otherwise the current solution is
// used. The secondary variables are used accordingly. This
// is required to compute the derivative of the storage term
// using the implicit Euler method.
const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[scvIdx];
storage[energyEqIdx] = 0;
// compute the energy storage
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
storage[energyEqIdx] +=
volVars.porosity()*(volVars.fluidState().density(phaseIdx) *
volVars.fluidState().internalEnergy(phaseIdx) *
volVars.fluidState().saturation(phaseIdx));
}
// heat capacity is already multiplied by the density
// of the porous material and the porosity in the problem file
storage[energyEqIdx] += volVars.temperature()*volVars.heatCapacity();
}
/*!
* \brief Evaluates the advective mass flux and the heat flux
* over a face of a sub-control volume and writes the result into
* the flux vector.
*
* \param flux The advective flux over the sub-control-volume face for each component
* \param fluxVars The flux variables at the current sub-control-volume face
*
* This method is called by compute flux (base class).
*/
void computeAdvectiveFlux(PrimaryVariables &flux,
const FluxVariables &fluxVars) const
{
// advective mass flux
ParentType::computeAdvectiveFlux(flux, fluxVars);
<