Commit 1567fb1c authored by Alexander Kissinger's avatar Alexander Kissinger
Browse files

- Adapted 3p3c to be used with the generic non-isothermal mode

- Deprecated the 3p3cni model
- Added effective thermal conductivity calculation in fluidmatrixinteractions for 3p
- Added thermal conductivity functions to the 3p fluidsystem h2oairmesitylen.hh and h2oairxylene.hh
- Upddated the 3p3cni test for the generic non-isothermal model
- Changed the reference vtus for the two tests, since results are slightly different due to different thermal conductivities


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@13427 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent dffb4df9
......@@ -32,6 +32,8 @@
#include <dumux/implicit/box/boxproperties.hh>
#include <dumux/implicit/cellcentered/ccproperties.hh>
#include <dumux/implicit/nonisothermal/niproperties.hh>
namespace Dumux
{
......@@ -47,6 +49,13 @@ NEW_TYPE_TAG(ThreePThreeC);
NEW_TYPE_TAG(BoxThreePThreeC, INHERITS_FROM(BoxModel, ThreePThreeC));
NEW_TYPE_TAG(CCThreePThreeC, INHERITS_FROM(CCModel, ThreePThreeC));
#ifndef DUMUX_3P3CNI_PROPERTIES_HH
//! The type tags for the corresponding non-isothermal problems
NEW_TYPE_TAG(ThreePThreeCNI, INHERITS_FROM(ThreePThreeC, NonIsothermal));
NEW_TYPE_TAG(BoxThreePThreeCNI, INHERITS_FROM(BoxModel, ThreePThreeCNI));
NEW_TYPE_TAG(CCThreePThreeCNI, INHERITS_FROM(CCModel, ThreePThreeCNI));
#endif
//////////////////////////////////////////////////////////////////
// Property tags
//////////////////////////////////////////////////////////////////
......
......@@ -38,9 +38,11 @@
#include "3p3cnewtoncontroller.hh"
#include "3p3clocalresidual.hh"
#include <dumux/implicit/nonisothermal/nipropertydefaults.hh>
#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
#include <dumux/material/spatialparams/implicitspatialparams.hh>
#include <dumux/material/fluidmatrixinteractions/3p/thermalconductivitysomerton3p.hh>
namespace Dumux
{
......@@ -148,6 +150,46 @@ SET_BOOL_PROP(ThreePThreeC, ProblemEnableGravity, true);
// (Nield, Bejan, Convection in porous media, 2006, p. 10)
SET_SCALAR_PROP(BoxModel, SpatialParamsForchCoeff, 0.55);
//! Somerton 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, Indices) Indices;
public:
typedef ThermalConductivitySomerton<Scalar, Indices> type;
};
//! temperature is already written by the isothermal model
SET_BOOL_PROP(ThreePThreeCNI, NiOutputLevel, 0);
//////////////////////////////////////////////////////////////////
// Property values for isothermal model required for the general non-isothermal model
//////////////////////////////////////////////////////////////////
// set isothermal Model
SET_TYPE_PROP(ThreePThreeCNI, IsothermalModel, ThreePThreeCModel<TypeTag>);
// set isothermal FluxVariables
SET_TYPE_PROP(ThreePThreeCNI, IsothermalFluxVariables, ThreePThreeCFluxVariables<TypeTag>);
//set isothermal VolumeVariables
SET_TYPE_PROP(ThreePThreeCNI, IsothermalVolumeVariables, ThreePThreeCVolumeVariables<TypeTag>);
//set isothermal LocalResidual
SET_TYPE_PROP(ThreePThreeCNI, IsothermalLocalResidual, ThreePThreeCLocalResidual<TypeTag>);
//set isothermal Indices
SET_PROP(ThreePThreeCNI, IsothermalIndices)
{
public:
typedef ThreePThreeCIndices<TypeTag, /*PVOffset=*/0> type;
};
//set isothermal NumEq
SET_INT_PROP(ThreePThreeCNI, IsothermalNumEq, 3);
}
}
......
......@@ -582,6 +582,12 @@ public:
// energy related quantities not contained in the fluid state
asImp_().updateEnergy_(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
// compute and set the enthalpy
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
Scalar h = Implementation::enthalpy_(fluidState_, paramCache, phaseIdx);
fluidState_.setEnthalpy(phaseIdx, h);
}
}
/*!
......@@ -720,6 +726,14 @@ protected:
bool isOldSol)
{ }
template<class ParameterCache>
static Scalar enthalpy_(const FluidState& fluidState,
const ParameterCache& paramCache,
const int phaseIdx)
{
return 0;
}
Scalar sw_, sg_, sn_, pg_, pw_, pn_;
Scalar moleFrac_[numPhases][numComponents];
......
......@@ -28,6 +28,8 @@
#ifndef DUMUX_3P3CNI_PROPERTIES_HH
#define DUMUX_3P3CNI_PROPERTIES_HH
#warning You are including the old 2pni model which will be removed after 2.6. See CHANGELOG for details.
#include <dumux/implicit/3p3c/3p3cproperties.hh>
namespace Dumux
......
// -*- 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 Relation for the saturation-dependent effective thermal conductivity
*/
#ifndef THERMALCONDUCTIVITY_SOMERTON_3P_HH
#define THERMALCONDUCTIVITY_SOMERTON_3P_HH
#include <algorithm>
namespace Dumux
{
struct SimpleThreePIndices
{
static const int wPhaseIdx = 0;
static const int nPhaseIdx = 1;
static const int gPhaseIdx = 2;
};
/*!
* \ingroup fluidmatrixinteractionslaws
*
* \brief Relation for the saturation-dependent effective thermal conductivity
*
* The Somerton method computes the thermal conductivity of dry and the wet soil material.
* It is extended here to a three phase system of water (w), NAPL (n) and gas (g).
* It uses a root function of the water and NAPL saturation to compute the
* effective thermal conductivity for a three-phase fluidsystem. The individual thermal
* conductivities are calculated as geometric mean of the thermal conductivity of the porous
* material and of the respective fluid phase.
*
* The material law is:
* \f[
\lambda_\text{eff} = \lambdaG_\text{eff} + \sqrt{(S_w)} \left(\lambdaW_\text{eff} - \lambdaG_\text{eff}\right) +
\sqrt{(S_n)} \left(\lambdaN_\text{eff} - \lambdaG_\text{eff}\right)
\f]
*
* with
* \f[
\lambdaW_\text{eff} = \lambda_{solid}^{\left(1-\phi\right)}*\lambda_w^\phi
\f]
* and
*
* \f[
\lambdaN_\text{eff} = \lambda_{solid}^{\left(1-\phi\right)}*\lambda_n^\phi.
\f]
*
* * \f[
\lambdaG_\text{eff} = \lambda_{solid}^{\left(1-\phi\right)}*\lambda_g^\phi.
\f]
*/
template<class Scalar, class Indices = SimpleThreePIndices>
class ThermalConductivitySomerton
{
public:
/*!
* \brief effective thermal conductivity \f$[W/(m K)]\f$ after Somerton (1974) extended for a three phase system
*
* \param volVars volume variables
* \param spatialParams spatial parameters
* \param element element (to be passed to spatialParams)
* \param fvGeometry fvGeometry (to be passed to spatialParams)
* \param scvIdx scvIdx (to be passed to spatialParams)
*
* \return effective thermal conductivity \f$[W/(m K)]\f$ after Somerton (1974)
*
* This gives an interpolation of the effective thermal conductivities of a porous medium
* filled with the water phase (w), a NAPL phase (n) and a gas phase (g).
* These two effective conductivities are computed as geometric mean of the solid and the
* fluid conductivities and interpolated with the square root of the wetting saturation.
*/
template<class VolumeVariables, class SpatialParams, class Element, class FVGeometry>
static Scalar effectiveThermalConductivity(const VolumeVariables& volVars,
const SpatialParams& spatialParams,
const Element& element,
const FVGeometry& fvGeometry,
int scvIdx)
{
Scalar sw = volVars.saturation(Indices::wPhaseIdx);
Scalar sn = volVars.saturation(Indices::nPhaseIdx);
Scalar lambdaW = volVars.thermalConductivityFluid(Indices::wPhaseIdx);
Scalar lambdaN = volVars.thermalConductivityFluid(Indices::nPhaseIdx);
Scalar lambdaG = volVars.thermalConductivityFluid(Indices::gPhaseIdx);
Scalar lambdaSolid = volVars.thermalConductivitySolid();
Scalar porosity = volVars.porosity();
return effectiveThermalConductivity(sw, sn, lambdaW, lambdaN, lambdaG, lambdaSolid, porosity);
}
/*!
* \brief effective thermal conductivity \f$[W/(m K)]\f$ after Somerton (1974)
*
* \param sw The saturation of the wetting phase
* \param sn The saturation of the wetting phase
* \param lambdaW the thermal conductivity of the water phase
* \param lambdaN the thermal conductivity of the NAPL phase
* \param lambdaG the thermal conductivity of the gas phase
* \param lambdaSolid the thermal conductivity of the solid phase
* \param porosity The porosity
*
* \return effective thermal conductivity \f$[W/(m K)]\f$ after Somerton (1974)
*/
static Scalar effectiveThermalConductivity(const Scalar sw,
const Scalar sn,
const Scalar lambdaW,
const Scalar lambdaN,
const Scalar lambdaG,
const Scalar lambdaSolid,
const Scalar porosity)
{
const Scalar satW = std::max<Scalar>(0.0, sw);
const Scalar satN = std::max<Scalar>(0.0, sn);
// const Scalar lSw = 1.8; //std::pow(lambdaSolid, (1.0 - porosity)) * std::pow(lambdaW, porosity);
// const Scalar lSn = 0.65; //std::pow(lambdaSolid, (1.0 - porosity)) * std::pow(lambdaN, porosity);
// const Scalar lSg = 0.35; //std::pow(lambdaSolid, (1.0 - porosity)) * std::pow(lambdaG, porosity);
// porosity weighted geometric mean
Scalar lSw = std::pow(lambdaSolid, (1.0 - porosity)) * std::pow(lambdaW, porosity);
Scalar lSn = std::pow(lambdaSolid, (1.0 - porosity)) * std::pow(lambdaN, porosity);
Scalar lSg = std::pow(lambdaSolid, (1.0 - porosity)) * std::pow(lambdaG, porosity);
Scalar lambdaEff = lSg + std::sqrt(satW) * (lSw - lSg) + std::sqrt(satN) * (lSn -lSg);
return lambdaEff;
}
};
}
#endif
......@@ -503,7 +503,33 @@ public:
static Scalar thermalConductivity(const FluidState &fluidState,
int phaseIdx)
{
DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OAirMesitylene::thermalConductivity()");
if (phaseIdx == wPhaseIdx){
const Scalar temperature = fluidState.temperature(phaseIdx) ;
const Scalar pressure = fluidState.pressure(phaseIdx);
return H2O::liquidThermalConductivity(temperature, pressure);
}
else if (phaseIdx == nPhaseIdx)
{
//Thermal Conductivity of p-Xylene
//taken from the Dortmund Data Bank
//see http://www.ddbst.de/en/EED/PCP/TCN_C176.php
const Scalar lambdaPureXylene = 0.13;
return lambdaPureXylene; // conductivity of p-xylene [W/(m K)]
}
else if (phaseIdx == gPhaseIdx)
{
//gPhaseIdx i.e. air
// Isobaric Properties for Nitrogen in: NIST Standard
// see http://webbook.nist.gov/chemistry/fluid/
// evaluated at p=.1 MPa, T=20°C
// Nitrogen: 0.025398
// Oxygen: 0.026105
// lambda_air is approximately 0.78*lambda_N2+0.22*lambda_O2
const Scalar lambdaPureAir = 0.0255535;
return lambdaPureAir; // conductivity of pure air [W/(m K)]
}
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
private:
......
......@@ -502,7 +502,33 @@ public:
static Scalar thermalConductivity(const FluidState &fluidState,
int phaseIdx)
{
DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OAirXylene::thermalConductivity()");
if (phaseIdx == wPhaseIdx){
const Scalar temperature = fluidState.temperature(phaseIdx) ;
const Scalar pressure = fluidState.pressure(phaseIdx);
return H2O::liquidThermalConductivity(temperature, pressure);
}
else if (phaseIdx == nPhaseIdx)
{
//Thermal Conductivity of p-Xylene
//taken from the Dortmund Data Bank
//see http://www.ddbst.de/en/EED/PCP/TCN_C176.php
const Scalar lambdaPureXylene = 0.13;
return lambdaPureXylene; // conductivity of p-xylene [W/(m K)]
}
else if (phaseIdx == gPhaseIdx)
{
//gPhaseIdx i.e. air
// Isobaric Properties for Nitrogen in: NIST Standard
// see http://webbook.nist.gov/chemistry/fluid/
// evaluated at p=.1 MPa, T=20°C
// Nitrogen: 0.025398
// Oxygen: 0.026105
// lambda_air is approximately 0.78*lambda_N2+0.22*lambda_O2
const Scalar lambdaPureAir = 0.0255535;
return lambdaPureAir; // conductivity of pure air [W/(m K)]
}
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
private:
......
......@@ -2,7 +2,7 @@ add_dumux_test(test_box3p3cni test_box3p3cni test_box3p3cni.cc
${CMAKE_SOURCE_DIR}/bin/runTest.sh
${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
${CMAKE_SOURCE_DIR}/test/references/column3p3cnibox-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/columnxylolbox-00060.vtu
${CMAKE_CURRENT_BINARY_DIR}/columnxylolbox-00058.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_box3p3cni)
add_dumux_test(test_cc3p3cni test_cc3p3cni test_cc3p3cni.cc
......
......@@ -31,7 +31,7 @@
#include <dumux/material/fluidsystems/h2oairxylenefluidsystem.hh>
#include <dumux/implicit/3p3cni/3p3cnimodel.hh>
#include <dumux/implicit/3p3c/3p3cmodel.hh>
#include <dumux/implicit/common/implicitporousmediaproblem.hh>
#include "columnxylolspatialparams.hh"
......
......@@ -116,6 +116,9 @@ public:
fineHeatCap_ = 850.;
coarseHeatCap_ = 84000.;
// heat conductivity of granite
lambdaSolid_ = 2.8;
// residual saturations
fineMaterialParams_.setSwr(0.12);
fineMaterialParams_.setSwrx(0.12);
......@@ -235,49 +238,21 @@ public:
* (1 - porosity(element, fvGeometry, scvIdx));
}
/*!
* \brief Calculate the heat flux \f$[W/m^2]\f$ through the
* rock matrix based on the temperature gradient \f$[K / m]\f$
* \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
*
* This is only required for non-isothermal models.
*
* \param heatFlux The resulting heat flux vector
* \param fluxVars The flux variables
* \param elemVolVars The volume variables
* \param tempGrad The temperature gradient
* \param element The current finite element
* \param fvGeometry The finite volume geometry of the current element
* \param faceIdx 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(DimVector &heatFlux,
const FluxVariables &fluxVars,
const ElementVolumeVariables &elemVolVars,
const DimVector &tempGrad,
const Element &element,
const FVElementGeometry &fvGeometry,
int faceIdx) const
Scalar thermalConductivitySolid(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
static const Scalar lDry = 0.35;
static const Scalar lsw1 = 1.8;
static const Scalar lsn1 = 0.65;
// arithmetic mean of the liquid saturation and the porosity
const int i = fluxVars.face().i;
const int j = fluxVars.face().j;
Scalar sw = std::max(0.0, (elemVolVars[i].saturation(wPhaseIdx) +
elemVolVars[j].saturation(wPhaseIdx)) / 2);
Scalar sn = std::max(0.0, (elemVolVars[i].saturation(nPhaseIdx) +
elemVolVars[j].saturation(nPhaseIdx)) / 2);
// the heat conductivity of the matrix. in general this is a
// tensorial value, but we assume isotropic heat conductivity.
Scalar heatCond = lDry + sqrt(sw) * (lsw1-lDry) + sqrt(sn) * (lsn1-lDry);
// the matrix heat flux is the negative temperature gradient
// times the heat conductivity.
heatFlux = tempGrad;
heatFlux *= -heatCond;
return lambdaSolid_;
}
private:
......@@ -299,6 +274,8 @@ private:
MaterialLawParams fineMaterialParams_;
MaterialLawParams coarseMaterialParams_;
Scalar lambdaSolid_;
};
}
......
......@@ -32,7 +32,7 @@
#include <dumux/material/fluidsystems/h2oairmesitylenefluidsystem.hh>
#include <dumux/implicit/3p3cni/3p3cnimodel.hh>
#include <dumux/implicit/3p3c/3p3cmodel.hh>
#include <dumux/implicit/common/implicitporousmediaproblem.hh>
#include "kuevettespatialparams.hh"
......
......@@ -114,6 +114,9 @@ public:
finePorosity_ = 0.42;
coarsePorosity_ = 0.42;
// heat conductivity of granite
lambdaSolid_ = 2.8;
// residual saturations
fineMaterialParams_.setSwr(0.12);
fineMaterialParams_.setSwrx(0.12);
......@@ -275,6 +278,21 @@ public:
heatFlux *= -heatCond;
}
/*!
* \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
*
* \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
*/
Scalar thermalConductivitySolid(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return lambdaSolid_;
}
private:
bool isFineMaterial_(const GlobalPosition &globalPos) const
{
......@@ -293,6 +311,8 @@ private:
MaterialLawParams fineMaterialParams_;
MaterialLawParams coarseMaterialParams_;
Scalar lambdaSolid_;
};
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment