Commit 766cb8d2 authored by Ned Coltman's avatar Ned Coltman
Browse files

[mm][fuelcell] Adapt to new material laws

parent 59925954
......@@ -56,3 +56,6 @@ SolidDensity = 1430 # [kg/m^3]
SolidThermalConductivity = 15.6 # [W/(m*K)]
SolidHeatCapacity = 710 # [J/(kg*K)]
[SpatialParams]
Swr = 0.05
Snr = 0.05
......@@ -28,12 +28,7 @@
#include <dumux/common/properties.hh>
#include <dumux/material/spatialparams/fv.hh>
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include "./material/regularizedacosta.hh"
#include "./material/acosta.hh"
#include "./material/thermalconductivityconstant.hh"
namespace Dumux {
......@@ -80,11 +75,8 @@ class FuelCellLectureSpatialParams: public FVSpatialParams<GetPropType<TypeTag,
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using EffectiveLaw = RegularizedAcosta<Scalar>;
using PcKrSwCurve = FluidMatrix::AcostaPcSwDefault<Scalar>;
public:
using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
using MaterialLawParams = typename MaterialLaw::Params;
using PermeabilityType = Scalar;
/*!
......@@ -98,23 +90,14 @@ public:
{
// intrinsic permeabilities
K_ = 5.2e-11; // Acosta absolute permeability diffusion layer
// porosities
porosity_ = 0.78; // Acosta porosity diffusion layer
// thermalconductivity
lambdaSolid_ = 15.6; // [W/(m*K)] Acosta thermal conductivity used in capillary pressure-saturation
// residual saturations
materialParams_.setSwr(0.05); // value for imbibition from Acosta
materialParams_.setSnr(0.05); // assumption
typename PcKrSwCurve::BasicParams params(-1168.75, 8.5, -0.2, -700, 0.0);
pcKrSwCurve_ = std::make_unique<PcKrSwCurve>(params);
// parameters for the Acosta saturation capillary pressure relation
materialParams_.setAcA(-1168.75); // imbition -1168.75; drainage -600;
materialParams_.setAcB(8.5); // imbition 8.5; drainage 25;
materialParams_.setAcC(-0.2); // imbition -0.2; drainage -16;
materialParams_.setAcD(-700); // imbition -700; drainage -3300;
materialParams_.setAcE(0.0); // imbition 0.0; drainage 800;
}
~FuelCellLectureSpatialParams()
......@@ -126,9 +109,7 @@ public:
* \param globalPos The global position
*/
PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
{
return K_;
}
{ return K_; }
/*!
* \brief Defines the porosity \f$[-]\f$ of the spatial parameters.
......@@ -154,19 +135,15 @@ public:
*/
template<class FluidSystem>
int wettingPhaseAtPos(const GlobalPosition& globalPos) const
{
return FluidSystem::phase0Idx;
}
{ return FluidSystem::phase0Idx; }
/*!
* \brief Returns the parameter object for the Brooks-Corey material law which depends on the position.
*
* \param globalPos The global position
*/
const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
{
return materialParams_;
}
auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const
{ return makeFluidMatrixInteraction(*pcKrSwCurve_); }
/*!
* \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
......@@ -181,9 +158,7 @@ public:
Scalar solidHeatCapacity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return 710; // specific heat capacity of diffusion layer Acosta [J / (kg K)]
}
{ return 710; }// specific heat capacity of diffusion layer Acosta [J / (kg K)]
/*!
* \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix.
......@@ -198,9 +173,7 @@ public:
Scalar solidDensity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return 1430; // density of ELAT [kg/m^3] Wöhr
}
{ return 1430; } // density of ELAT [kg/m^3] Wöhr
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid.
......@@ -214,15 +187,13 @@ public:
Scalar solidThermalConductivity(const Element &element,
const FVElementGeometry &fvGeometry,
const SubControlVolume& scv) const
{
return lambdaSolid_;
}
{ return lambdaSolid_; }
private:
Scalar K_;
Scalar porosity_;
Scalar eps_ = 1e-6;
MaterialLawParams materialParams_;
std::unique_ptr<const PcKrSwCurve> pcKrSwCurve_;
Scalar lambdaSolid_;
};
......
......@@ -25,17 +25,15 @@
#ifndef ACOSTA_HH
#define ACOSTA_HH
#include "acostaparams.hh"
#include <algorithm>
#include <cmath>
#include <cassert>
#include <vector>
#include <dumux/common/exceptions.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/spline.hh>
#include <dumux/common/optionalscalar.hh>
#include <dumux/material/fluidmatrixinteractions/2p/materiallaw.hh>
namespace Dumux
{
namespace Dumux::FluidMatrix {
/*!
*
......@@ -48,14 +46,62 @@ namespace Dumux
*
* \see AcostaParams
*/
template <class ScalarT, class ParamsT = AcostaPcSwParams<ScalarT> >
class AcostaPcSw
{
public:
template<class Scalar>
struct Params
{
Params(Scalar acA, Scalar acB, Scalar acC, Scalar acD, Scalar acE)
: acA_(acA)
, acB_(acB)
, acC_(acC)
, acD_(acD)
, acE_(acE)
{}
Scalar acA() const { return acA_; }
void setAcA(Scalar a) { acA_ = a; }
Scalar acB() const { return acB_; }
void setAcB(Scalar b) { acB_ = b; }
Scalar acC() const { return acC_; }
void setAcC(Scalar c) { acC_ = c; }
Scalar acD() const { return acD_; }
void setAcD(Scalar d) { acD_ = d; }
typedef ParamsT Params;
typedef typename Params::Scalar Scalar;
Scalar acE() const { return acE_; }
void setAcE(Scalar e) { acE_ = e; }
bool operator== (const Params& p) const
{
return Dune::FloatCmp::eq(acA(), p.acA(), 1e-6)
&& Dune::FloatCmp::eq(acB(), p.acB(), 1e-6)
&& Dune::FloatCmp::eq(acC(), p.acC(), 1e-6)
&& Dune::FloatCmp::eq(acD(), p.acD(), 1e-6)
&& Dune::FloatCmp::eq(acE(), p.acE(), 1e-6);
}
private:
Scalar acA_, acB_, acC_, acD_, acE_;
};
/*!
* \brief Construct from a subgroup from the global parameter tree
* \note This will give you nice error messages if a mandatory parameter is missing
*/
template<class Scalar = double>
static Params<Scalar> makeParams(const std::string& paramGroup)
{
const auto acA = getParamFromGroup<Scalar>(paramGroup, "AcA");
const auto acB = getParamFromGroup<Scalar>(paramGroup, "AcB");
const auto acC = getParamFromGroup<Scalar>(paramGroup, "AcC");
const auto acD = getParamFromGroup<Scalar>(paramGroup, "AcD");
const auto acE = getParamFromGroup<Scalar>(paramGroup, "AcE");
return {acA, acB, acC, acD, acE};
}
/*!
* \brief The capillary pressure-saturation curve according to Acosta.
*
......@@ -71,7 +117,8 @@ public:
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*/
static Scalar pc(const Params &params, Scalar swe)
template<class Scalar>
static Scalar pc(Scalar swe, const Params<Scalar>& params)
{
assert(0 <= swe && swe <= 1);
Scalar pc = exp(params.acB() * swe + params.acC());
......@@ -94,9 +141,10 @@ public:
* is constructed accordingly. Afterwards the values are set there, too.
* \return The effective saturation of the wetting phase
*/
static Scalar sw(const Params &params, Scalar pc)
template<class Scalar>
static Scalar sw(Scalar pc, const Params<Scalar>& params)
{
DUNE_THROW(Dune::NotImplemented, "Acosta::sw(params, pc)");
DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::sw(params, pc)");
}
/*!
......@@ -111,7 +159,8 @@ public:
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*/
static Scalar dpc_dsw(const Params &params, Scalar swe)
template<class Scalar>
static Scalar dpc_dsw(Scalar swe, const Params<Scalar>& params)
{
assert(0 <= swe && swe <= 1);
......@@ -134,9 +183,10 @@ public:
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*/
static Scalar dsw_dpc(const Params &params, Scalar pc)
template<class Scalar>
static Scalar dsw_dpc(Scalar pc, const Params<Scalar>& params)
{
DUNE_THROW(Dune::NotImplemented, "Acosta::dsw_dpc(params, pc)");
DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dsw_dpc(params, pc)");
}
/*!
......@@ -148,7 +198,8 @@ public:
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too. */
static Scalar krw(const Params &params, Scalar swe)
template<class Scalar>
static Scalar krw(Scalar swe, const Params<Scalar>& params)
{
assert(0 <= swe && swe <= 1);
Scalar krw;
......@@ -180,9 +231,10 @@ public:
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*/
static Scalar dkrw_dsw(const Params &params, Scalar swe)
template<class Scalar>
static Scalar dkrw_dsw(Scalar swe, const Params<Scalar>& params)
{
DUNE_THROW(Dune::NotImplemented, "Acosta::dkrw_dsw(params, swe)");
DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dkrw_dsw(params, swe)");
};
/*!
......@@ -195,7 +247,8 @@ public:
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*/
static Scalar krn(const Params &params, Scalar swe)
template<class Scalar>
static Scalar krn(Scalar swe, const Params<Scalar>& params)
{
assert(0 <= swe && swe <= 1);
......@@ -229,13 +282,14 @@ public:
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
*/
static Scalar dkrn_dsw(const Params &params, Scalar swe)
template<class Scalar>
static Scalar dkrn_dsw(Scalar swe, const Params<Scalar>& params)
{
DUNE_THROW(Dune::NotImplemented, "Acosta::dkrn_dsw(params, swe)");
DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dkrn_dsw(params, swe)");
}
private:
template<class Scalar>
static Scalar evaluateIntegral_(const std::vector<double>& intVector, const Scalar swe)
{
Scalar satPos = (integrationSteps_ - 1) * swe;
......@@ -257,23 +311,226 @@ private:
static std::vector<double> intDraSw1_;
};
template <class ScalarT, class ParamsT>
std::vector<double> AcostaPcSw<ScalarT, ParamsT>::intImb0Sw_ = {0.0000000000e+00, 1.4448654771e-08, 2.2819077599e-08, 2.7262869143e-08, 2.9465634517e-08, 3.0504079857e-08,
template <class Scalar>
class AcostaRegularization
{
public:
template<class S>
struct Params
{
Params(S thresholdSw = 1e-3)
: thresholdSw_(thresholdSw)
{}
S thresholdSw() const
{ return thresholdSw_; }
bool operator== (const Params& p) const
{ return Dune::FloatCmp::eq(thresholdSw(), p.thresholdSw(), 1e-6);}
private:
S thresholdSw_;
};
template<class MaterialLaw, class BaseParams, class EffToAbsParams>
void init(const MaterialLaw* m, const BaseParams& bp, const EffToAbsParams& etap, const Params<Scalar>& p)
{
thresholdSw_ = p.thresholdSw();
initParameters_(m, bp, thresholdSw_);
}
/*!
* \brief A regularized Acosta capillary pressure-saturation
* curve.
*
* regularized part:
* - low saturation: extend the \f$p_c(S_w)\f$ curve with the slope at the regularization point (i.e. no kink).
* - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$ by a straight line (yes, there is a kink :-( ).
*
* For the non-regularized part:
*
* \copydetails AcostaPcSw::pc()
*/
OptionalScalar<Scalar> pc(Scalar swe) const
{
// make sure that the capilary pressure observes a derivative != 0 for 'illegal' saturations. This is
// required for example by newton solvers (if the derivative is calculated numerically) in order to get the
// saturation moving to the right direction if it temporarily is in an 'illegal' range.
if (acE_ == 0.0) //imbibition
{
if (swe < thresholdSw_)
return pcsweLow_ + mLow_*( swe - thresholdSw_);
else if (swe >= (1.0 - thresholdSw_))
return pcsweHigh_ + mHigh_*(swe - (1.0 - thresholdSw_));
}
else //drainage
{
if (swe <= thresholdSw_)
return pcsweLow_ + mLow_*(swe - thresholdSw_);
else if (swe >= (1.0 - thresholdSw_))
return pcsweHigh_ + mHigh_*(swe - (1.0 - thresholdSw_));
}
// if the effective saturation is in an 'reasonable'
// range, we use the real Acosta law...
return {};
}
/*!
* \brief A regularized Acosta saturation-capillary pressure curve.
*
* function does not exist yet!
*/
OptionalScalar<Scalar> sw(Scalar pc) const
{
DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::sw(bp, pc)");
}
/*!
* \brief A regularized version of the partial derivative
* of the \f$p_c(\overline S_w)\f$ w.r.t. effective saturation
* according to Acosta.
*
* regularized part:
* - low saturation: use the slope of the regularization point (i.e. no kink).
* - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$
* by a straight line and use that slope (yes, there is a kink :-( ).
*
* For the non-regularized part:
*
* \copydetails AcostaPcSw::dpc_dsw()
*/
OptionalScalar<Scalar> dpc_dsw(Scalar swe) const
{
// derivative of the regualarization
if (swe < 0)
return m0_;
else if (swe >= (1.0 - thresholdSw_))
return mHigh_;
return {};
}
/*!
* \brief The partial derivative of the effective
* saturation to the capillary pressure according to Acosta.
*
* function does not exist yet!
*/
OptionalScalar<Scalar> dsw_dpc(Scalar pc) const
{ DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dsw_dpc(pc, bp)"); }
/*!
* \brief Regularized version of the relative permeability
* for the wetting phase of
* the medium implied by the Acosta
* parameterization.
*
* regularized part:
* - below \f$ \overline S_w =0\f$: set relative permeability to zero
* - above \f$ \overline S_w =1\f$: set relative permeability to one
* - between \f$ 0.95 \leq \overline S_w \leq 1\f$: use a spline as interpolation
*
* For not-regularized part: \copydetails AcostaPcSw::krw()
*/
OptionalScalar<Scalar> krw(Scalar swe) const
{
if (swe <= 0.0)
return 0.0;
else if (swe >= 1.0)
return 1.0;
return {};
}
/*!
* \brief The derivative of the relative permeability for the
* wetting phase in regard to the wetting saturation of the
* medium implied by the Acosta parameterization.
*/
OptionalScalar<Scalar> dkrw_dsw(Scalar swe) const
{ DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dkrw_dsw(swe, bp)"); };
/*!
* \brief Regularized version of the relative permeability
* for the nonwetting phase of
* the medium implied by the Acosta
* parameterization.
*
* regularized part:
* - below \f$ \overline S_w =0\f$: set relative permeability to zero
* - above \f$ \overline S_w =1\f$: set relative permeability to one
* - for \f$ 0 \leq \overline S_w \leq 0.05 \f$: use a spline as interpolation
* \copydetails AcostaPcSw::krn()
*/
OptionalScalar<Scalar> krn(Scalar swe) const
{
if (swe >= 1.0)
return 0.0;
else if (swe <= 0.0)
return 1.0;
return {};
}
/*!
* \brief The derivative of the relative permeability for the
* nonwetting phase in regard to the wetting saturation of
* the medium as implied by the Acosta
* parameterization.
*/
OptionalScalar<Scalar> dkrn_dsw(Scalar swe) const
{ DUNE_THROW(Dune::NotImplemented, "AcostaPcSw::dkrn_dsw(swe, bp)"); };
private:
template<class MaterialLaw, class BaseParams>
void initParameters_(const MaterialLaw* m, const BaseParams& bp, const Scalar thresholdSw)
{
thresholdSw_ = thresholdSw;
m0_ = AcostaPcSw::dpc_dsw(0.0, bp);
mLow_ = AcostaPcSw::dpc_dsw(thresholdSw_, bp);
pcsweLow_ = AcostaPcSw::pc(thresholdSw_, bp);
mHigh_ = AcostaPcSw::dpc_dsw((1.0-thresholdSw_), bp);
pcsweHigh_ = AcostaPcSw::pc((1.0-thresholdSw_), bp);
acE_ = bp.acE();
}
Scalar thresholdSw_;
Scalar m0_;
Scalar mLow_;
Scalar pcsweLow_;
Scalar mHigh_;
Scalar pcsweHigh_;
Scalar acE_;
};
/*!
* \ingroup Fluidmatrixinteractions
* \brief A default configuration for using the VanGenuchten material law
*/
template<typename Scalar = double>
using AcostaPcSwDefault = TwoPMaterialLaw<Scalar, AcostaPcSw, AcostaRegularization<Scalar>, TwoPEffToAbsDefaultPolicy>;
std::vector<double> AcostaPcSw::intImb0Sw_ = {0.0000000000e+00, 1.4448654771e-08, 2.2819077599e-08, 2.7262869143e-08, 2.9465634517e-08, 3.0504079857e-08,
3.0976893747e-08, 3.1187237078e-08, 3.1279418104e-08, 3.1319431548e-08, 3.1336696665e-08, 3.1344118639e-08, 3.1347301932e-08, 3.1348665344e-08,
3.1349248803e-08, 3.1349498362e-08, 3.1349605071e-08, 3.1349650691e-08, 3.1349670192e-08, 3.1349678527e-08, 3.1349682090e-08};
template <class ScalarT, class ParamsT>
std::vector<double> AcostaPcSw<ScalarT, ParamsT>::intImbSw1_ = {3.1349682090e-08, 1.6901027319e-08, 8.5306044914e-09, 4.0868129471e-09, 1.8840475735e-09, 8.4560223338e-10,
std::vector<double> AcostaPcSw::intImbSw1_ = {3.1349682090e-08, 1.6901027319e-08, 8.5306044914e-09, 4.0868129471e-09, 1.8840475735e-09, 8.4560223338e-10,
3.7278834340e-10, 1.6244501261e-10, 7.0263986151e-11, 3.0250542321e-11, 1.2985425390e-11, 5.5634511134e-12, 2.3801587239e-12, 1.0167462066e-12,
4.3328717089e-13, 1.8372876392e-13, 7.7019533374e-14, 3.1399697336e-14, 1.1898564029e-14, 3.5629091876e-15, 0.0000000000e+00};
template <class ScalarT, class ParamsT>
std::vector<double> AcostaPcSw<ScalarT, ParamsT>::intDra0Sw_ = {0.0000000000e+00, 9.0286303155e-11, 1.0440347151e-09, 5.3755677032e-09, 2.0922261595e-08, 7.4676307329e-08,
std::vector<double> AcostaPcSw::intDra0Sw_ = {0.0000000000e+00, 9.0286303155e-11, 1.0440347151e-09, 5.3755677032e-09, 2.0922261595e-08, 7.4676307329e-08,
2.8041525597e-07, 1.3662599439e-06, 2.4661875855e-05, 9.2186966645e-01, 9.2188585863e-01, 9.2189413795e-01, 9.2189728035e-01, 9.2189773354e-01,
9.2189777056e-01, 9.2189777338e-01, 9.2189777360e-01, 9.2189777360e-01, 9.2189777359e-01, 9.2189777356e-01, 9.2189777362e-01};
template <class ScalarT, class ParamsT>
std::vector<double> AcostaPcSw<ScalarT, ParamsT>::intDraSw1_ = {9.2189777362e-01, 9.2189777350e-01, 9.2189777259e-01, 9.2189776823e-01, 9.2189775271e-01, 9.2189769898e-01,
std::vector<double> AcostaPcSw::intDraSw1_ = {9.2189777362e-01, 9.2189777350e-01, 9.2189777259e-01, 9.2189776823e-01, 9.2189775271e-01, 9.2189769898e-01,
9.2189749322e-01, 9.2189640736e-01, 9.2187311174e-01, 2.8107160428e-05, 1.1914958378e-05, 3.6356677731e-06, 4.9325835742e-07, 4.0074681584e-08,
3.0250536841e-09, 2.3593717317e-10, 1.8918448614e-11, 1.5374841055e-12, 1.2497552223e-13, 9.4684316652e-15, 0.0000000000e+00};
}
......
// -*- 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 Specification of the material parameters
* for the Acosta constitutive relations.
*/
#ifndef ACOSTA_PARAMS_HH
#define ACOSTA_PARAMS_HH
namespace Dumux
{
/*!
*
* \brief Specification of the material parameters
* for the Acosta constitutive relations.
*/
template<class ScalarT>
class AcostaPcSwParams
{
public:
typedef ScalarT Scalar;
AcostaPcSwParams()
{}
AcostaPcSwParams(Scalar acA, Scalar acB, Scalar acC, Scalar acD, Scalar acE)
{ setAcA(acA); setAcB(acB); setAcC(acC); setAcD(acD); setAcE(acE); };
/*!
* \ Acosta's curve.
*/
Scalar acA() const
{ return acA_; }
/*!
* \ Acosta's curve.
*/
void setAcA(Scalar a)
{ acA_ = a; }
/*!