diff --git a/dumux/implicit/co2ni/co2nivolumevariables.hh b/dumux/implicit/co2ni/co2nivolumevariables.hh
index f4a3894548884be094543478457ffff1d2ea813d..70af07fb855d3fae040f378b8a4d74a0fb2e99f1 100644
--- a/dumux/implicit/co2ni/co2nivolumevariables.hh
+++ b/dumux/implicit/co2ni/co2nivolumevariables.hh
@@ -87,6 +87,14 @@ public:
Scalar heatCapacity() const
{ return heatCapacity_; };
+ /*!
+ * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of the fluid phase in
+ * the sub-control volume.
+ */
+ Scalar thermalConductivity(const int phaseIdx) const
+ { return FluidSystem::thermalConductivity(this->fluidState_, phaseIdx); };
+
+
protected:
// this method gets called by the parent class. since this method
// is protected, we are friends with our parent..
diff --git a/test/implicit/co2/heterogeneousproblem.hh b/test/implicit/co2/heterogeneousproblem.hh
index 95504a05ea7a98ad1a357fc4b46d750f7ea27888..a1634854cd1594331e4c117d5dcb3d8822d74ded 100644
--- a/test/implicit/co2/heterogeneousproblem.hh
+++ b/test/implicit/co2/heterogeneousproblem.hh
@@ -106,7 +106,7 @@ SET_BOOL_PROP(HeterogeneousProblem, VtkAddVelocity, false);
/*!
* \ingroup CO2Model
- * \ingroup BoxTestProblems
+ * \ingroup ImplicitTestProblems
* \brief Problem where CO2 is injected under a low permeable layer in a depth of 1200m.
*
* The domain is sized 200m times 100m and consists of four layers, a
@@ -122,7 +122,7 @@ SET_BOOL_PROP(HeterogeneousProblem, VtkAddVelocity, false);
* These boundary ids can be imported into the problem where the boundary conditions can then be assigned accordingly.
*
* To run the simulation execute the following line in shell:
- * <tt>./test_co2 </tt>
+ * <tt>./test_ccco2 </tt> or <tt>./test_boxco2 </tt>
*/
template <class TypeTag >
class HeterogeneousProblem : public ImplicitPorousMediaProblem<TypeTag>
diff --git a/test/implicit/co2/heterogeneousspatialparameters.hh b/test/implicit/co2/heterogeneousspatialparameters.hh
index 7a60d6ef50f2f31e68e8addb1498f2dbb9f1ec12..81202442bab03faced8a54d1cf47530cbfb4374a 100644
--- a/test/implicit/co2/heterogeneousspatialparameters.hh
+++ b/test/implicit/co2/heterogeneousspatialparameters.hh
@@ -24,7 +24,7 @@
* \file
*
* \brief Definition of the spatial parameters for the injection
- * problem which uses the isothermal CO2 box model
+ * problem which uses the non-isothermal or isothermal CO2 box or cc model
*/
#ifndef DUMUX_HETEROGENEOUS_SPATIAL_PARAMS_HH
@@ -70,7 +70,7 @@ public:
* \ingroup CO2Model
* \ingroup BoxTestProblems
* \brief Definition of the spatial parameters for the injection
- * problem which uses the isothermal CO2 box model
+ * problem which uses the non-isothermal or isothermal CO2NI box or cc model
*/
template<class TypeTag>
class HeterogeneousSpatialParams : public ImplicitSpatialParams<TypeTag>
@@ -114,14 +114,17 @@ public:
{
/*
* Layer Index Setup:
- * barrierTop = 4
- * barrierMiddle = 5
- * reservoir = 6
+ * barrierTop = 1
+ * barrierMiddle = 2
+ * reservoir = 3
*/
barrierTop_ = 1;
barrierMiddle_ = 2;
reservoir_ = 3;
+ // heat conductivity of granite
+ lambdaSolid_ = 2.8;
+
//Set the permeability tensor for the layers
//Scalar anisotropy = 0.1;
//for (int i = 0; i < dim; i++)
@@ -165,7 +168,6 @@ public:
int elemIdx = (*gridPtr_)->leafView().indexSet().index(*elemIt);
int param = (*gridPtr_).parameters(*elemIt)[0];
paramIdx_[elemIdx] = param;
- //std::cout<<"param: "<<paramIdx_[elemIdx]<<std::endl;
}
}
@@ -251,52 +253,22 @@ public:
}
/*!
- * \brief Calculate the heat flux \f$[W/m^2]\f$ through the
- * rock matrix based on the temperature gradient \f$[K / m]\f$
- *
- * This is only required for non-isothermal models.
+ * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
*
- * \param heatFlux The resulting heat flux vector
- * \param fluxDat The flux variables
- * \param vDat The volume variables
- * \param tempGrad The temperature gradient
- * \param element The current finite element
- * \param fvElemGeom The finite volume geometry of the current element
- * \param scvfIdx The local index of the sub-control volume face where
- * the matrix heat flux should be calculated
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the heat capacity needs to be defined
*/
- void matrixHeatFlux(Vector &heatFlux,
- const FluxVariables &fluxDat,
- const ElementVolumeVariables &vDat,
- const Vector &tempGrad,
- const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvfIdx) const
+ Scalar thermalConductivitySolid(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
{
- static const Scalar lWater = 0.6;
- static const Scalar lGranite = 2.8;
-
- // arithmetic mean of the liquid saturation and the porosity
- const int i = fvElemGeom.subContVolFace[scvfIdx].i;
- const int j = fvElemGeom.subContVolFace[scvfIdx].j;
- Scalar Sl = std::max<Scalar>(0.0, (vDat[i].saturation(lPhaseIdx) +
- vDat[j].saturation(lPhaseIdx)) / 2);
- Scalar poro = (porosity(element, fvElemGeom, i) +
- porosity(element, fvElemGeom, j)) / 2;
-
- Scalar lsat = pow(lGranite, (1-poro)) * pow(lWater, poro);
- Scalar ldry = pow(lGranite, (1-poro));
-
- // the heat conductivity of the matrix. in general this is a
- // tensorial value, but we assume isotropic heat conductivity.
- Scalar heatCond = ldry + sqrt(Sl) * (ldry - lsat);
-
- // the matrix heat flux is the negative temperature gradient
- // times the heat conductivity.
- heatFlux = tempGrad;
- heatFlux *= -heatCond;
+ return lambdaSolid_;
}
+
+
private:
@@ -312,6 +284,7 @@ private:
Scalar barrierTopK_;
Scalar barrierMiddleK_;
Scalar reservoirK_;
+ Scalar lambdaSolid_;
MaterialLawParams materialParams_;
diff --git a/test/implicit/co2ni/heterogeneousproblemni.hh b/test/implicit/co2ni/heterogeneousproblemni.hh
index 1376d85d0d10f2c67d2dcacb4758ec6ac90ef32c..d03d7bfa6a32b795dbafd94536ad2437b49a4c5d 100644
--- a/test/implicit/co2ni/heterogeneousproblemni.hh
+++ b/test/implicit/co2ni/heterogeneousproblemni.hh
@@ -36,9 +36,8 @@
#include <dumux/material/fluidsystems/brineco2fluidsystem.hh>
#include <dumux/implicit/common/implicitporousmediaproblem.hh>
#include <dumux/implicit/box/intersectiontovertexbc.hh>
+#include <test/implicit/co2/heterogeneousspatialparameters.hh>
-
-#include "heterogeneousspatialparametersni.hh"
#include "heterogeneousco2tables.hh"
#define ISOTHERMAL 0
diff --git a/test/implicit/co2ni/heterogeneousspatialparametersni.hh b/test/implicit/co2ni/heterogeneousspatialparametersni.hh
deleted file mode 100644
index 3a21d6ae5e5c3cea2233589e04328f058b06925b..0000000000000000000000000000000000000000
--- a/test/implicit/co2ni/heterogeneousspatialparametersni.hh
+++ /dev/null
@@ -1,323 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2008-2009 by Klaus Mosthaf *
- * Copyright (C) 2008-2009 by Andreas Lauser *
- * Institute for Modelling Hydraulic and Environmental Systems *
- * University of Stuttgart, Germany *
- * email: <givenname>.<name>@iws.uni-stuttgart.de *
- * *
- * 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 Definition of the spatial parameters for the injection
- * problem which uses the non-isothermal CO2 box model
- */
-
-#ifndef DUMUX_HETEROGENEOUS_NI_SPATIAL_PARAMS_HH
-#define DUMUX_HETEROGENEOUS_NI_SPATIAL_PARAMS_HH
-
-#include <dumux/material/spatialparams/implicitspatialparams.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
-
-#include <dumux/implicit/co2/co2model.hh>
-
-namespace Dumux
-{
-
-//forward declaration
-template<class TypeTag>
-class HeterogeneousSpatialParams;
-
-namespace Properties
-{
-// The spatial parameters TypeTag
-NEW_TYPE_TAG(HeterogeneousSpatialParams);
-
-// Set the spatial parameters
-SET_TYPE_PROP(HeterogeneousSpatialParams, SpatialParams, Dumux::HeterogeneousSpatialParams<TypeTag>);
-
-// Set the material Law
-SET_PROP(HeterogeneousSpatialParams, MaterialLaw)
-{
-private:
- // define the material law which is parameterized by effective
- // saturations
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw;
-public:
- // define the material law parameterized by absolute saturations
- typedef EffToAbsLaw<EffMaterialLaw> type;
-};
-}
-
-/*!
- * \ingroup CO2NIModel
- * \ingroup BoxTestProblems
- * \brief Definition of the spatial parameters for the injection
- * problem which uses the non-isothermal CO2NI box model
- */
-template<class TypeTag>
-class HeterogeneousSpatialParams : public ImplicitSpatialParams<TypeTag>
-{
- typedef ImplicitSpatialParams<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef Dune::GridPtr<Grid> GridPointer;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename Grid::ctype CoordScalar;
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-
- enum {
- dim=GridView::dimension,
- dimWorld=GridView::dimensionworld,
-
- lPhaseIdx = FluidSystem::lPhaseIdx
- };
-
- typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
- typedef Dune::FieldVector<CoordScalar,dimWorld> Vector;
-
- typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
-
-public:
- typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
- typedef typename MaterialLaw::Params MaterialLawParams;
-
- /*!
- * \brief The constructor
- *
- * \param gridView The grid view
- */
- HeterogeneousSpatialParams(const GridView &gridView)
- : ParentType(gridView)
- {
- /*
- * Layer Index Setup:
- * barrierTop = 1
- * barrierMiddle = 2
- * reservoir = 3
- */
- barrierTop_ = 1;
- barrierMiddle_ = 2;
- reservoir_ = 3;
-
- //Set the permeability tensor for the layers
- //Scalar anisotropy = 0.1;
- //for (int i = 0; i < dim; i++)
- barrierTopK_ = 1e-17; //sqm
- //barrierTopK_[dim-1][dim-1] = barrierTopK_[0][0]*anisotropy;
-
- //for (int i = 0; i < dim; i++)
- barrierMiddleK_ = 1e-15; //sqm
- //barrierMiddleK_[dim-1][dim-1] = barrierMiddleK_[0][0]*anisotropy;
-
- // for (int i = 0; i < dim; i++)
- reservoirK_ = 1e-14; //sqm
- //reservoirK_[dim-1][dim-1] = reservoirK_[0][0]*anisotropy;
-
- //Set the effective porosity of the layers
- barrierTopPorosity_ = 0.001;
- barrierMiddlePorosity_ = 0.05;
- reservoirPorosity_ = 0.2;
-
-
- // Same material parameters for every layer
- materialParams_.setSwr(0.2);
- materialParams_.setSwr(0.05);
- materialParams_.setLambda(2.0);
- materialParams_.setPe(1e4);
- }
-
- ~HeterogeneousSpatialParams()
- {}
-
- void setParams(GridPointer *gridPtr)
- {
- gridPtr_ = gridPtr;
- int numElems = (*gridPtr_)->leafView().size(0);
- paramIdx_.resize(numElems);
-
- ElementIterator elemIt = (*gridPtr_)->leafView().template begin<0>();
- const ElementIterator elemEndIt = (*gridPtr_)->leafView().template end<0>();
- for (; elemIt != elemEndIt; ++elemIt)
- {
- int elemIdx = (*gridPtr_)->leafView().indexSet().index(*elemIt);
- int param = (*gridPtr_).parameters(*elemIt)[0];
- paramIdx_[elemIdx] = param;
- }
-
- }
-
- /*!
- * \brief Apply the intrinsic permeability tensor to a pressure
- * potential gradient.
- *
- * \param element The current finite element
- * \param fvElemGeom The current finite volume geometry of the element
- * \param scvIdx The index of the sub-control volume
- */
- const Scalar intrinsicPermeability(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
- {
- int elemIdx = (*gridPtr_)->leafView().indexSet().index(element); //Get the global index of the element
-
- if (paramIdx_[elemIdx] == barrierTop_)
- return barrierTopK_;
- else if (paramIdx_[elemIdx] == barrierMiddle_)
- return barrierMiddleK_;
- else
- return reservoirK_;
- }
-
- /*!
- * \brief Define the porosity \f$[-]\f$ of the spatial parameters
- *
- * \param element The finite element
- * \param fvElemGeom The finite volume geometry
- * \param scvIdx The local index of the sub-control volume where
- * the porosity needs to be defined
- */
- Scalar porosity(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
- {
- int elemIdx = (*gridPtr_)->leafView().indexSet().index(element); //Get the global index of the element
-
- if (paramIdx_[elemIdx] == barrierTop_)
- return barrierTopPorosity_;
- else if (paramIdx_[elemIdx] == barrierMiddle_)
- return barrierMiddlePorosity_;
- else
- return reservoirPorosity_;
- }
-
-
- /*!
- * \brief return the parameter object for the Brooks-Corey material law which depends on the position
- *
- * \param element The current finite element
- * \param fvElemGeom The current finite volume geometry of the element
- * \param scvIdx The index of the sub-control volume
- */
- const MaterialLawParams& materialLawParams(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
- {
-
- return materialParams_;
- }
-
- /*!
- * \brief Returns the heat capacity \f$[J/m^3 K]\f$ of the rock matrix.
- *
- * This is only required for non-isothermal models.
- *
- * \param element The finite element
- * \param fvElemGeom The finite volume geometry
- * \param scvIdx The local index of the sub-control volume where
- * the heat capacity needs to be defined
- */
- double heatCapacity(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
- {
- return
- 790 // specific heat capacity of granite [J / (kg K)]
- * 2700 // density of granite [kg/m^3]
- * (1 - porosity(element, fvElemGeom, scvIdx));
- }
-
- /*!
- * \brief Calculate the heat flux \f$[W/m^2]\f$ through the
- * rock matrix based on the temperature gradient \f$[K / m]\f$
- *
- * This is only required for non-isothermal models.
- *
- * \param heatFlux The resulting heat flux vector
- * \param fluxDat The flux variables
- * \param vDat The volume variables
- * \param tempGrad The temperature gradient
- * \param element The current finite element
- * \param fvElemGeom The finite volume geometry of the current element
- * \param scvfIdx The local index of the sub-control volume face where
- * the matrix heat flux should be calculated
- */
- void matrixHeatFlux(Vector &heatFlux,
- const FluxVariables &fluxDat,
- const ElementVolumeVariables &vDat,
- const Vector &tempGrad,
- const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvfIdx) const
- {
- static const Scalar lWater = 0.6;
- static const Scalar lGranite = 2.8;
-
- // arithmetic mean of the liquid saturation and the porosity
- const int i = fvElemGeom.subContVolFace[scvfIdx].i;
- const int j = fvElemGeom.subContVolFace[scvfIdx].j;
- Scalar Sl = std::max<Scalar>(0.0, (vDat[i].saturation(lPhaseIdx) +
- vDat[j].saturation(lPhaseIdx)) / 2);
- Scalar poro = (porosity(element, fvElemGeom, i) +
- porosity(element, fvElemGeom, j)) / 2;
-
- Scalar lsat = pow(lGranite, (1-poro)) * pow(lWater, poro);
- Scalar ldry = pow(lGranite, (1-poro));
-
- // the heat conductivity of the matrix. in general this is a
- // tensorial value, but we assume isotropic heat conductivity.
- Scalar heatCond = ldry + sqrt(Sl) * (ldry - lsat);
-
- // the matrix heat flux is the negative temperature gradient
- // times the heat conductivity.
- heatFlux = tempGrad;
- heatFlux *= -heatCond;
- }
-
-private:
-
-
- int barrierTop_;
- int barrierMiddle_;
- int reservoir_;
-
-
- Scalar barrierTopPorosity_;
- Scalar barrierMiddlePorosity_;
- Scalar reservoirPorosity_;
-
- Scalar barrierTopK_;
- Scalar barrierMiddleK_;
- Scalar reservoirK_;
-
- MaterialLawParams materialParams_;
-
- GridPointer *gridPtr_;
- std::vector<int> paramIdx_;
-};
-
-}
-
-#endif