Commit 8465d28a authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[nonisothermal] improve naming and meaning of solid thermal parameters

This implements FS#216.

The "heatCapacity" function of the spatial parameters and the volume
variables for the implicit nonisothermal models was a misnomer, since it
returned an effective quantity, namely, 
heatCapacity*density*(1 - porosity) in [J/(K m^3)].
Except for mpnc, which resulted in an additional inconsistency.

Corresponding to the decision documented in FS#216, this patch renames
the function to "solidHeatCapacity" and returns always the "true"
(non-effective) heat capacity in [J/(kg K)]. This requires an additional
function "solidDensity" which returns the mass density of the porous
matrix. Moreover, the functions "thermalConductivitySolid/Fluid" are
renamed to "solid/fluidThermalConductivity". The decision to prepend
with "solid/fluid" rather than to append is motivated by consistency
with components and fluid systems, where "gas" and "liquid" are always
prepended to the corresponding function names.

_Beware_: this change breaks compatibility. You have to adapt your
spatial parameters such that they offer functions "solidHeatCapacity",
"solidDensity" and "solidThermalConductivity".

Reviewed by Alex.



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14070 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 09d3f6b1
......@@ -429,8 +429,8 @@ class BoxFVElementGeometry
{0, 2, 0, 1, 3, 2}
};
static const int edgeToFacePyramid[2][8] = {
{1, 2, 3, 4, 1, 3, 4, 2},
{0, 0, 0, 0, 3, 2, 1, 4}
{0, 2, 3, 0, 1, 3, 4, 2},
{1, 0, 0, 4, 3, 2, 1, 4}
};
static const int edgeToFacePrism[2][9] = {
{1, 0, 2, 0, 3, 2, 4, 1, 4},
......
......@@ -19,7 +19,7 @@
/*!
* \file
*
* \brief Adaption of the BOX scheme to the two-phase two-component flow model.
* \brief velocity output for implicit (porous media) models
*/
#ifndef DUMUX_IMPLICIT_VELOCITYOUTPUT_HH
#define DUMUX_IMPLICIT_VELOCITYOUTPUT_HH
......@@ -47,6 +47,7 @@ class ImplicitVelocityOutput
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::ctype CoordScalar;
......@@ -63,6 +64,7 @@ class ImplicitVelocityOutput
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef typename Dune::ReferenceElements<CoordScalar, dim> ReferenceElements;
typedef typename Dune::ReferenceElement<CoordScalar, dim> ReferenceElement;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
......@@ -78,41 +80,31 @@ public:
{
// check, if velocity output can be used (works only for cubes so far)
velocityOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity);
if (velocityOutput_ && isBox)
{
cellNum_.assign(problem_.gridView().size(dofCodim), 0);
}
ElementIterator eIt = problem_.gridView().template begin<0>();
ElementIterator eEndIt = problem_.gridView().template end<0>();
for (; eIt != eEndIt; ++eIt)
if (velocityOutput_)
{
if (eIt->geometry().type().isCube() == false)
{
velocityOutput_ = false;
}
if (velocityOutput_ && isBox)
if (isBox)
{
FVElementGeometry fvGeometry;
fvGeometry.update(problem_.gridView(), *eIt);
cellNum_.assign(problem_.gridView().size(dofCodim), 0);
// transform vertex velocities from local to global coordinates
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
ElementIterator eIt = problem_.gridView().template begin<0>();
ElementIterator eEndIt = problem_.gridView().template end<0>();
for (; eIt != eEndIt; ++eIt)
{
FVElementGeometry fvGeometry;
fvGeometry.update(problem_.gridView(), *eIt);
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
{
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
int vIdxGlobal = problem_.vertexMapper().subIndex(*eIt, scvIdx, dofCodim);
int vIdxGlobal = problem_.vertexMapper().subIndex(*eIt, scvIdx, dofCodim);
#else
int vIdxGlobal = problem_.vertexMapper().map(*eIt, scvIdx, dofCodim);
int vIdxGlobal = problem_.vertexMapper().map(*eIt, scvIdx, dofCodim);
#endif
cellNum_[vIdxGlobal] += 1;
cellNum_[vIdxGlobal] += 1;
}
}
}
}
if (velocityOutput_ != GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity))
std::cout << "ATTENTION: Velocity output only works for cubes and is set to false for simplices\n";
}
bool enableOutput()
......@@ -124,16 +116,17 @@ public:
void calculateVelocity(VelocityVector& velocity,
const ElementVolumeVariables& elemVolVars,
const FVElementGeometry& fvGeometry,
const Element& element,
const Element& element,
int phaseIdx)
{
if (velocityOutput_)
{
// calculate vertex velocities
GlobalPosition tmpVelocity(0.0);
Dune::GeometryType geomType = element.geometry().type();
const ReferenceElement &referenceElement
= ReferenceElements::general(geomType);
const Dune::FieldVector<Scalar, dim>& localPos =
ReferenceElements::general(element.geometry().type()).position(0, 0);
const Dune::FieldVector<Scalar, dim>& localPos
= referenceElement.position(0, 0);
// get the transposed Jacobian of the element mapping
const typename Element::Geometry::JacobianTransposed jacobianT2 =
......@@ -141,8 +134,8 @@ public:
if (isBox)
{
typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > SCVVelocities;
SCVVelocities scvVelocities(pow(2,dim));
typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > ScvVelocities;
ScvVelocities scvVelocities(fvGeometry.numScv);
scvVelocities = 0;
for (int fIdx = 0; fIdx < fvGeometry.numScvf; fIdx++)
......@@ -164,17 +157,17 @@ public:
GlobalPosition localNormal(0);
jacobianT1.mv(globalNormal, localNormal);
// note only works for cubes
const Scalar localArea = pow(2,-(dim-1));
localNormal /= localNormal.two_norm();
// Get the Darcy velocities. The Darcy velocities are divided by the area of the subcontrolvolume
// face in the reference element.
// area of the subcontrolvolume face in the reference element
Scalar localArea = scvfReferenceArea_(geomType, fIdx);
// get the volume flux divided by the area of the
// subcontrolvolume face in the reference element
Scalar flux = fluxVars.volumeFlux(phaseIdx) / localArea;
// transform the normal Darcy velocity into a vector
tmpVelocity = localNormal;
// transform the volume flux into a velocity vector
GlobalPosition tmpVelocity = localNormal;
tmpVelocity *= flux;
scvVelocities[fluxVars.face().i] += tmpVelocity;
......@@ -200,7 +193,11 @@ public:
}
else
{
Dune::FieldVector<Scalar, 2*dim> scvVelocities(0.0);
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
std::vector<Scalar> scvfFluxes(element.subEntities(1), 0);
#else
std::vector<Scalar> scvfFluxes(element.template count<1>(), 0);
#endif
int fIdxInner = 0;
IntersectionIterator isEndIt = problem_.gridView().iend(element);
......@@ -217,8 +214,7 @@ public:
fIdxInner,
elemVolVars);
Scalar flux = fluxVars.volumeFlux(phaseIdx);
scvVelocities[fIdx] += flux;
scvfFluxes[fIdx] += fluxVars.volumeFlux(phaseIdx);
fIdxInner++;
}
......@@ -230,14 +226,68 @@ public:
fIdx,
elemVolVars,true);
Scalar flux = fluxVars.volumeFlux(phaseIdx);
scvVelocities[fIdx] = flux;
scvfFluxes[fIdx] = fluxVars.volumeFlux(phaseIdx);
}
}
// Correct boundary fluxes in case of Neumann conditions.
// They are simply set to an average of the inner fluxes. In
// this general setting, it would be very difficult to
// calculate correct phase, i.e., volume, fluxes from arbitrary
// Neumann conditions.
if (element.hasBoundaryIntersections())
{
IntersectionIterator isEndIt = problem_.gridView().iend(element);
for (IntersectionIterator isIt = problem_.gridView().ibegin(element);
isIt != isEndIt; ++isIt)
{
if (isIt->boundary())
{
BoundaryTypes bcTypes;
problem_.boundaryTypes(bcTypes, *isIt);
if (bcTypes.hasNeumann())
{
int fIdx = isIt->indexInInside();
// cubes
if (dim == 1 || geomType.isCube()){
int fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
}
// simplices
else if (geomType.isSimplex()) {
scvfFluxes[fIdx] = 0;
}
}
}
}
}
Dune::FieldVector < Scalar, dim > refVelocity(0);
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (scvVelocities[2*i + 1] - scvVelocities[2*i]);
Dune::FieldVector < Scalar, dim > refVelocity;
// cubes
if (dim == 1 || geomType.isCube()){
for (int i = 0; i < dim; i++)
{
refVelocity[i] = 0.5 * (scvfFluxes[2*i + 1] - scvfFluxes[2*i]);
}
}
// simplices
else if (geomType.isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -scvfFluxes[dim - 1 - dimIdx];
for (int fIdx = 0; fIdx < dim + 1; fIdx++)
{
refVelocity[dimIdx] += scvfFluxes[fIdx]/(dim + 1);
}
}
}
// 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented,
"velocity output for cell-centered and prism/pyramid");
}
Dune::FieldVector<Scalar, dimWorld> scvVelocity(0);
jacobianT2.mtv(refVelocity, scvVelocity);
......@@ -255,6 +305,63 @@ public:
} // velocity output
}
private:
// The area of a subcontrolvolume face in a reference element.
// The 3d non-cube values have been calculated with quadrilateralArea3D
// of boxfvelementgeometry.hh.
static Scalar scvfReferenceArea_(Dune::GeometryType geomType, int fIdx)
{
if (dim == 1 || geomType.isCube())
{
return 1.0/(1 << (dim-1));
}
else if (geomType.isTriangle())
{
static const Scalar faceToArea[] = {0.372677996249965,
0.372677996249965,
0.235702260395516};
return faceToArea[fIdx];
}
else if (geomType.isTetrahedron())
{
static const Scalar faceToArea[] = {0.102062072615966,
0.102062072615966,
0.058925565098879,
0.102062072615966,
0.058925565098879,
0.058925565098879};
return faceToArea[fIdx];
}
else if (geomType.isPyramid())
{
static const Scalar faceToArea[] = {0.130437298687488,
0.130437298687488,
0.130437298687488,
0.130437298687488,
0.150923085635624,
0.1092906420717,
0.1092906420717,
0.0781735959970571};
return faceToArea[fIdx];
}
else if (geomType.isPrism())
{
static const Scalar faceToArea[] = {0.166666666666667,
0.166666666666667,
0.166666666666667,
0.186338998124982,
0.186338998124982,
0.117851130197758,
0.186338998124982,
0.186338998124982,
0.117851130197758};
return faceToArea[fIdx];
}
else {
DUNE_THROW(Dune::NotImplemented, "scvf area for unknown GeometryType");
}
}
protected:
const Problem& problem_;
bool velocityOutput_;
......
......@@ -195,13 +195,13 @@ protected:
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].fluidState().saturation(wPhaseIdx),
elemVolVars[i].thermalConductivity(wPhaseIdx),
elemVolVars[i].thermalConductivity(nPhaseIdx),
problem.spatialParams().thermalConductivitySolid(element, fvGeometry, i),
problem.spatialParams().solidThermalConductivity(element, fvGeometry, i),
problem.spatialParams().porosity(element, fvGeometry, i));
lambdaJ =
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].fluidState().saturation(wPhaseIdx),
elemVolVars[j].thermalConductivity(wPhaseIdx),
elemVolVars[j].thermalConductivity(nPhaseIdx),
problem.spatialParams().thermalConductivitySolid(element, fvGeometry, j),
problem.spatialParams().solidThermalConductivity(element, fvGeometry, j),
problem.spatialParams().porosity(element, fvGeometry, j));
}
else
......@@ -214,7 +214,7 @@ protected:
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].fluidState().saturation(wPhaseIdx),
elemVolVars[i].thermalConductivity(wPhaseIdx),
elemVolVars[i].thermalConductivity(nPhaseIdx),
problem.spatialParams().thermalConductivitySolid(elementI, fvGeometryI, 0),
problem.spatialParams().solidThermalConductivity(elementI, fvGeometryI, 0),
problem.spatialParams().porosity(elementI, fvGeometryI, 0));
const Element & elementJ = *fvGeometry.neighbors[j];
......@@ -225,7 +225,7 @@ protected:
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].fluidState().saturation(wPhaseIdx),
elemVolVars[j].thermalConductivity(wPhaseIdx),
elemVolVars[j].thermalConductivity(nPhaseIdx),
problem.spatialParams().thermalConductivitySolid(elementJ, fvGeometryJ, 0),
problem.spatialParams().solidThermalConductivity(elementJ, fvGeometryJ, 0),
problem.spatialParams().porosity(elementJ, fvGeometryJ, 0));
}
......
......@@ -248,13 +248,13 @@ protected:
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].fluidState().saturation(wPhaseIdx),
elemVolVars[i].thermalConductivity(wPhaseIdx),
elemVolVars[i].thermalConductivity(nPhaseIdx),
problem.spatialParams().thermalConductivitySolid(element, fvGeometry, i),
problem.spatialParams().solidThermalConductivity(element, fvGeometry, i),
problem.spatialParams().porosity(element, fvGeometry, i));
lambdaJ =
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].fluidState().saturation(wPhaseIdx),
elemVolVars[j].thermalConductivity(wPhaseIdx),
elemVolVars[j].thermalConductivity(nPhaseIdx),
problem.spatialParams().thermalConductivitySolid(element, fvGeometry, j),
problem.spatialParams().solidThermalConductivity(element, fvGeometry, j),
problem.spatialParams().porosity(element, fvGeometry, j));
}
else
......@@ -267,7 +267,7 @@ protected:
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].fluidState().saturation(wPhaseIdx),
elemVolVars[i].thermalConductivity(wPhaseIdx),
elemVolVars[i].thermalConductivity(nPhaseIdx),
problem.spatialParams().thermalConductivitySolid(elementI, fvGeometryI, 0),
problem.spatialParams().solidThermalConductivity(elementI, fvGeometryI, 0),
problem.spatialParams().porosity(elementI, fvGeometryI, 0));
const Element & elementJ = *fvGeometry.neighbors[j];
......@@ -278,7 +278,7 @@ protected:
ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].fluidState().saturation(wPhaseIdx),
elemVolVars[j].thermalConductivity(wPhaseIdx),
elemVolVars[j].thermalConductivity(nPhaseIdx),
problem.spatialParams().thermalConductivitySolid(elementJ, fvGeometryJ, 0),
problem.spatialParams().solidThermalConductivity(elementJ, fvGeometryJ, 0),
problem.spatialParams().porosity(elementJ, fvGeometryJ, 0));
}
......
......@@ -162,9 +162,9 @@ public:
// heat stored in the rock matrix
storage[energyEqIdx] +=
volVars.fluidState().temperature(/*phaseIdx=*/0)
* volVars.densitySolid()
* volVars.solidDensity()
* (1.0 - volVars.porosity())
* volVars.heatCapacity();
* volVars.solidHeatCapacity();
}
/*!
* \brief Calculate the storage for all mass balance equations
......
......@@ -121,9 +121,9 @@ public:
else if(phaseIdx == sPhaseIdx) {
// heat stored in the rock matrix
storage[energyEq0Idx+phaseIdx] += volVars.temperature(phaseIdx) *
volVars.densitySolid() *
volVars.solidDensity() *
(1.-volVars.porosity()) *
volVars.heatCapacity();
volVars.solidHeatCapacity();
}
else
DUNE_THROW(Dune::NotImplemented,
......@@ -463,9 +463,9 @@ public:
// heat stored in the rock matrix
storage[energyEqSolidIdx] +=
volVars.temperature(temperatureSolidIdx)
* volVars.densitySolid()
* volVars.solidDensity()
* (1.0 - volVars.porosity())
* volVars.heatCapacity();
* volVars.solidHeatCapacity();
if (!std::isfinite(storage[energyEqSolidIdx]))
DUNE_THROW(NumericalProblem, "Calculated non-finite storage");
......
......@@ -190,13 +190,13 @@ public:
Valgrind::SetUndefined(*this);
// heat capacities of the fluids plus the porous medium
heatCapacity_ =
problem.spatialParams().heatCapacity(element, fvGeometry, scvIdx);
Valgrind::CheckDefined(heatCapacity_);
solidHeatCapacity_ =
problem.spatialParams().solidHeatCapacity(element, fvGeometry, scvIdx);
Valgrind::CheckDefined(solidHeatCapacity_);
densitySolid_ =
problem.spatialParams().densitySolid(element, fvGeometry, scvIdx);
Valgrind::CheckDefined(densitySolid_);
solidDensity_ =
problem.spatialParams().solidDensity(element, fvGeometry, scvIdx);
Valgrind::CheckDefined(solidDensity_);
// set the enthalpies
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
......@@ -207,11 +207,15 @@ public:
}
/*!
* \brief Returns the total heat capacity [J/(K m^3)] of the rock matrix in
* \brief Returns the total heat capacity [J/(kg K)] of the rock matrix in
* the sub-control volume.
*/
Scalar solidHeatCapacity() const
{ return solidHeatCapacity_; }
Scalar heatCapacity() const
{ return heatCapacity_; }
DUNE_DEPRECATED_MSG("use solidHeatCapacity() instead")
{ return solidHeatCapacity(); }
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of the fluid phase in
......@@ -224,9 +228,12 @@ public:
* \brief Returns the total density of the given solid phase [kg / m^3] in
* the sub-control volume.
*/
Scalar densitySolid() const
{ return densitySolid_; }
Scalar solidDensity() const
{ return solidDensity_; }
Scalar densitySolid() const
DUNE_DEPRECATED_MSG("use solidDensity() instead")
{ return solidDensity(); }
/*!
* \brief If running under valgrind this produces an error message
......@@ -234,13 +241,13 @@ public:
*/
void checkDefined() const
{
Valgrind::CheckDefined(heatCapacity_);
Valgrind::CheckDefined(densitySolid_);
Valgrind::CheckDefined(solidHeatCapacity_);
Valgrind::CheckDefined(solidDensity_);
}
protected:
Scalar heatCapacity_;
Scalar densitySolid_;
Scalar solidHeatCapacity_;
Scalar solidDensity_;
Scalar thermalConductivity_[numPhases] ;
};
......
......@@ -135,24 +135,24 @@ public:
const unsigned int scvIdx,
const Problem & problem)
{
heatCapacity_ =
problem.spatialParams().heatCapacity(element, fvGeometry, scvIdx);
Valgrind::CheckDefined(heatCapacity_);
solidHeatCapacity_ =
problem.spatialParams().solidHeatCapacity(element, fvGeometry, scvIdx);
Valgrind::CheckDefined(solidHeatCapacity_);
for(int phaseIdx =0; phaseIdx<numPhases; ++phaseIdx){
thermalConductivityFluid_[phaseIdx] =
fluidThermalConductivity_[phaseIdx] =
FluidSystem::thermalConductivity(fluidState, paramCache, phaseIdx);
}
Valgrind::CheckDefined(thermalConductivityFluid_);
Valgrind::CheckDefined(fluidThermalConductivity_);
densitySolid_ =
problem.spatialParams().densitySolid(element, fvGeometry, scvIdx);
Valgrind::CheckDefined(densitySolid_);
solidDensity_ =
problem.spatialParams().solidDensity(element, fvGeometry, scvIdx);
Valgrind::CheckDefined(solidDensity_);
thermalConductivitySolid_ =
problem.spatialParams().thermalConductivitySolid(element, fvGeometry, scvIdx);
Valgrind::CheckDefined(thermalConductivitySolid_);
solidThermalConductivity_ =
problem.spatialParams().solidThermalConductivity(element, fvGeometry, scvIdx);
Valgrind::CheckDefined(solidThermalConductivity_);
// set the enthalpies
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
......@@ -163,11 +163,15 @@ public:
}
/*!
* \brief Returns the total heat capacity [J/(K m^3)] of the rock matrix in
* \brief Returns the total heat capacity [J/(kg K)] of the rock matrix in
* the sub-control volume.
*/
Scalar solidHeatCapacity() const
{ return solidHeatCapacity_; }
Scalar heatCapacity() const
{ return heatCapacity_; }
DUNE_DEPRECATED_MSG("use solidHeatCapacity() instead")
{ return solidHeatCapacity(); }
/*!
* \brief Returns the temperature in fluid / solid phase(s)
......@@ -181,18 +185,26 @@ public:
* \brief Returns the total density of the given solid phase [kg / m^3] in
* the sub-control volume.
*/
Scalar solidDensity() const
{ return solidDensity_; }
Scalar densitySolid() const
{ return densitySolid_; }
DUNE_DEPRECATED_MSG("use solidDensity() instead")
{ return solidDensity(); }
/*!
* \brief Returns the conductivity of the given solid phase [kg / m^3] in
* \brief Returns the conductivity of the given solid phase [W/(m K)] in
* the sub-control volume.
*/
Scalar solidThermalConductivity() const
{ return solidThermalConductivity_; }
Scalar thermalConductivitySolid() const
{ return thermalConductivitySolid_; }
DUNE_DEPRECATED_MSG("use solidThermalConductivity() instead")
{ return solidThermalConductivity(); }
/*!
* \brief Returns the conductivity of the given fluid [kg / m^3] in
* \brief Returns the conductivity of the given fluid [W//m K)] in
* the sub-control volume.
*
* \param phaseIdx The local index of the phases
......@@ -200,9 +212,9 @@ public:
Scalar thermalConductivity(const unsigned int phaseIdx) const
{
if(phaseIdx == wPhaseIdx or phaseIdx == nPhaseIdx )
return thermalConductivityFluid_[phaseIdx];
return fluidThermalConductivity_[phaseIdx];
else if (phaseIdx == sPhaseIdx )
return thermalConductivitySolid_;
return solidThermalConductivity_;
else
DUNE_THROW(Dune::NotImplemented,