mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
big cleanup
timemanager API, restart, boxmodel are affected
This commit is contained in:
Andreas Lauser
committed by
Andreas Lauser
Andreas Lauser
parent
c07de5d462
commit
63ebb939bf
@@ -1,4 +1,4 @@
|
||||
// $Id: fullermethod.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: h2o_n2.hh 3783 2010-06-24 11:33:53Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: henryiapws.hh 3783 2010-06-24 11:33:53Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: component.hh 3824 2010-07-13 13:30:02Z lauser $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser
|
||||
* Institute of Hydraulic Engineering *
|
||||
@@ -21,8 +21,6 @@
|
||||
#ifndef DUMUX_COMPONENT_HH
|
||||
#define DUMUX_COMPONENT_HH
|
||||
|
||||
#include <dune/common/exceptions.hh>
|
||||
|
||||
namespace Dumux
|
||||
{
|
||||
|
||||
@@ -33,6 +31,10 @@ template <class Scalar, class Implementation>
|
||||
class Component
|
||||
{
|
||||
public:
|
||||
static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
|
||||
Scalar pressMin, Scalar pressMax, unsigned nPress)
|
||||
{ Dune::dwarn << "No init routine defined - make shure that this is not necessary!" << std::endl; }
|
||||
|
||||
/*!
|
||||
* \brief A human readable name for the compoent.
|
||||
*/
|
||||
|
||||
@@ -106,7 +106,7 @@ public:
|
||||
/*!
|
||||
* \brief The density [kg/m^3] of H2 at a given pressure and temperature.
|
||||
*/
|
||||
static Scalar density(Scalar temperature, Scalar pressure)
|
||||
static Scalar gasDensity(Scalar temperature, Scalar pressure)
|
||||
{
|
||||
// Assume an ideal gas
|
||||
return IdealGas::density(molarMass(), temperature, pressure);
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: h2o.hh 3783 2010-06-24 11:33:53Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: common.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009-2010 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: region1.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009-2010 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: region2.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009-2010 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: region4.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009-2010 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: n2.hh 3783 2010-06-24 11:33:53Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -116,7 +116,7 @@ public:
|
||||
*
|
||||
* \todo density liquid oxygen
|
||||
*/
|
||||
static Scalar density(Scalar temperature, Scalar pressure)
|
||||
static Scalar gasDensity(Scalar temperature, Scalar pressure)
|
||||
{
|
||||
// Assume an ideal gas
|
||||
return IdealGas::density(molarMass(), temperature, pressure);
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: simplednapl.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: simpleh2o.hh 3783 2010-06-24 11:33:53Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: tabulatedcomponent.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: unit.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: brookscorey.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2008 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: brookscoreyparams.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: efftoabslaw.hh 3811 2010-07-05 12:59:17Z pnuske $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: efftoabslawparams.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: linearmaterial.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: linearmaterialparams.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: regularizedbrookscorey.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: regularizedbrookscoreyparams.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: regularizedvangenuchten.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: regularizedvangenuchtenparams.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: vangenuchten.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: vangenuchtenparams.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: fluidstate.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser
|
||||
* Institute of Hydraulic Engineering *
|
||||
@@ -117,8 +117,8 @@ public:
|
||||
*
|
||||
* Unit: [Pa] = [N/m^2]
|
||||
*/
|
||||
Scalar partialPressure(int componentIdx) const
|
||||
{ DUNE_THROW(Dune::NotImplemented, "FluidState::partialPressure()"); }
|
||||
Scalar fugacity(int componentIdx) const
|
||||
{ DUNE_THROW(Dune::NotImplemented, "FluidState::fugacity()"); }
|
||||
|
||||
/*!
|
||||
* \brief Return the total pressure of the gas phase.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: 2p_system.hh 3784 2010-06-24 13:43:57Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: h2o_n2_system.hh 3824 2010-07-13 13:30:02Z lauser $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser
|
||||
* Institute of Hydraulic Engineering *
|
||||
@@ -91,57 +91,6 @@ public:
|
||||
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Given the gas phase's composition, temperature and
|
||||
* pressure, compute the partial presures of all components
|
||||
* in [Pa] and assign it to the FluidState.
|
||||
*
|
||||
* This is required for models which cannot calculate the the
|
||||
* partial pressures of the components in the gas phase from the
|
||||
* degasPressure(). To use this method, the FluidState must have a
|
||||
* setPartialPressure(componentIndex, pressure) method.
|
||||
*/
|
||||
template <class FluidState>
|
||||
static void computePartialPressures(Scalar temperature,
|
||||
Scalar pg,
|
||||
FluidState &fluidState)
|
||||
{
|
||||
Scalar X1 = fluidState.massFrac(gPhaseIdx, H2OIdx);
|
||||
|
||||
// We use the newton method for this. For the initial value we
|
||||
// assume all components to be an ideal gas
|
||||
Scalar pH2O = fluidState.moleFrac(gPhaseIdx, H2OIdx) * pg;
|
||||
Scalar eps = pg*1e-9;
|
||||
|
||||
Scalar deltaP = pH2O;
|
||||
Valgrind::CheckDefined(pH2O);
|
||||
Valgrind::CheckDefined(deltaP);
|
||||
for (int i = 0; i < 5 && std::abs(deltaP/pg) > 1e-9; ++i) {
|
||||
Scalar f =
|
||||
H2O::gasDensity(temperature, pH2O)*(X1 - 1) +
|
||||
X1*N2::gasDensity(temperature, pg - pH2O);
|
||||
|
||||
Scalar df_dp;
|
||||
df_dp =
|
||||
H2O::gasDensity(temperature, pH2O + eps)*(X1 - 1) +
|
||||
X1*N2::gasDensity(temperature, pg - (pH2O + eps));
|
||||
df_dp -=
|
||||
H2O::gasDensity(temperature, pH2O - eps)*(X1 - 1) +
|
||||
X1*N2::gasDensity(temperature, pg - (pH2O - eps));
|
||||
df_dp /=
|
||||
2*eps;
|
||||
|
||||
deltaP = - f/df_dp;
|
||||
|
||||
pH2O += deltaP;
|
||||
Valgrind::CheckDefined(pH2O);
|
||||
Valgrind::CheckDefined(deltaP);
|
||||
}
|
||||
|
||||
fluidState.setPartialPressure(H2OIdx, pH2O);
|
||||
fluidState.setPartialPressure(N2Idx, pg - pH2O);
|
||||
};
|
||||
|
||||
/*!
|
||||
* \brief Given a phase's composition, temperature, pressure, and
|
||||
* the partial pressures of all components, return its
|
||||
@@ -153,16 +102,26 @@ public:
|
||||
Scalar pressure,
|
||||
const FluidState &fluidState)
|
||||
{
|
||||
if (phaseIdx == lPhaseIdx)
|
||||
return liquidPhaseDensity_(temperature,
|
||||
pressure,
|
||||
fluidState.moleFrac(lPhaseIdx, H2OIdx),
|
||||
fluidState.moleFrac(lPhaseIdx, N2Idx));
|
||||
else if (phaseIdx == gPhaseIdx)
|
||||
return gasPhaseDensity_(temperature,
|
||||
pressure,
|
||||
fluidState.moleFrac(gPhaseIdx, H2OIdx),
|
||||
fluidState.moleFrac(gPhaseIdx, N2Idx));
|
||||
if (phaseIdx == lPhaseIdx) {
|
||||
// See: Ochs 2008
|
||||
// \todo: proper citation
|
||||
Scalar rholH2O = H2O::liquidDensity(temperature, pressure);
|
||||
Scalar clH2O = rholH2O/H2O::molarMass();
|
||||
|
||||
// this assumes each nitrogen molecule displaces exactly one
|
||||
// water molecule in the liquid
|
||||
return
|
||||
clH2O*(H2O::molarMass()*fluidState.moleFrac(lPhaseIdx, H2OIdx)
|
||||
+
|
||||
N2::molarMass()*fluidState.moleFrac(lPhaseIdx, N2Idx));
|
||||
}
|
||||
else if (phaseIdx == gPhaseIdx) {
|
||||
Scalar fugH2O = fluidState.fugacity(H2OIdx);
|
||||
Scalar fugN2 = fluidState.fugacity(N2Idx);
|
||||
return
|
||||
H2O::gasDensity(temperature, fugH2O) +
|
||||
N2::gasDensity(temperature, fugN2);
|
||||
}
|
||||
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
|
||||
}
|
||||
|
||||
@@ -337,23 +296,46 @@ public:
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the fugacity coefficient of a component in a
|
||||
* \brief Returns the activity coefficient of a component in a
|
||||
* phase.
|
||||
*
|
||||
* For solutions with only traces in a solvent this boils down to
|
||||
* the inverse Henry constant for the solutes and the partial
|
||||
* pressure for the solvent.
|
||||
* We define the activity coefficent \f$\gamma_{\alpha,\kappa}\f$
|
||||
* of component \f$\kappa\f$ by the following equation:
|
||||
* \f[ f_\kappa = p_\alpha \gamma_{\alpha,\kappa} \f]
|
||||
* where \f$f_\kappa\f$ is the component's fugacity and \f$p_\alpha\f$
|
||||
* is the phase' pressure
|
||||
*
|
||||
* For liquids with very low miscibility this boils down to the
|
||||
* inverse Henry constant for the solutes and the partial pressure
|
||||
* for the solvent.
|
||||
*
|
||||
* For ideal gases this is equivalent to the gas pressure, in real
|
||||
* gases it is the gas pressure times the component's fugacity
|
||||
* coefficient.
|
||||
*/
|
||||
template <class FluidState>
|
||||
static Scalar fugacityCoeff(int phaseIdx,
|
||||
static Scalar activityCoeff(int phaseIdx,
|
||||
int compIdx,
|
||||
Scalar temperature,
|
||||
Scalar pressure,
|
||||
const FluidState &state)
|
||||
{
|
||||
if (phaseIdx != lPhaseIdx)
|
||||
DUNE_THROW(Dune::NotImplemented,
|
||||
"Fugacities in the gas phase");
|
||||
if (phaseIdx == gPhaseIdx) {
|
||||
return pressure;
|
||||
Scalar fugH2O = std::max(1e-3, state.fugacity(H2OIdx));
|
||||
Scalar fugN2 = std::max(1e-3, state.fugacity(N2Idx));
|
||||
Scalar cH2O = H2O::gasDensity(temperature, fugH2O) / H2O::molarMass();
|
||||
Scalar cN2 = N2::gasDensity(temperature, fugN2) / N2::molarMass();
|
||||
|
||||
Scalar alpha = (fugH2O + fugN2)/pressure;
|
||||
|
||||
if (compIdx == H2OIdx)
|
||||
return fugH2O/(alpha*cH2O/(cH2O + cN2));
|
||||
else if (compIdx == N2Idx)
|
||||
return fugN2/(alpha*cN2/(cH2O + cN2));
|
||||
|
||||
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
|
||||
}
|
||||
|
||||
switch (compIdx) {
|
||||
case H2OIdx: return H2O::vaporPressure(temperature);
|
||||
@@ -362,59 +344,6 @@ public:
|
||||
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Given a component's pressure and temperature, return its
|
||||
* density in a phase [kg/m^3].
|
||||
*/
|
||||
static Scalar componentDensity(int phaseIdx,
|
||||
int compIdx,
|
||||
Scalar temperature,
|
||||
Scalar pressure)
|
||||
{
|
||||
if (phaseIdx == lPhaseIdx) {
|
||||
if (compIdx == H2OIdx)
|
||||
return H2O::liquidDensity(temperature, pressure);
|
||||
else if (compIdx == N2Idx)
|
||||
return N2::liquidDensity(temperature, pressure);
|
||||
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
|
||||
}
|
||||
else if (phaseIdx == gPhaseIdx) {
|
||||
if (compIdx == H2OIdx) {
|
||||
return H2O::gasDensity(temperature, pressure);
|
||||
}
|
||||
else if (compIdx == N2Idx)
|
||||
return N2::gasDensity(temperature, pressure);
|
||||
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
|
||||
}
|
||||
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
|
||||
};
|
||||
|
||||
/*!
|
||||
* \brief Given a component's density and temperature, return the
|
||||
* corresponding pressure in a phase [Pa].
|
||||
*/
|
||||
static Scalar componentPressure(int phaseIdx,
|
||||
int compIdx,
|
||||
Scalar temperature,
|
||||
Scalar density)
|
||||
{
|
||||
if (phaseIdx == lPhaseIdx) {
|
||||
if (compIdx == H2OIdx)
|
||||
return H2O::liquidPressure(temperature, density);
|
||||
else if (compIdx == N2Idx)
|
||||
return N2::liquidPressure(temperature, density);
|
||||
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
|
||||
}
|
||||
else if (phaseIdx == gPhaseIdx) {
|
||||
if (compIdx == H2OIdx)
|
||||
return H2O::gasPressure(temperature, density);
|
||||
else if (compIdx == N2Idx)
|
||||
return N2::gasPressure(temperature, density);
|
||||
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
|
||||
}
|
||||
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
|
||||
};
|
||||
|
||||
/*!
|
||||
* \brief Given a phase's composition, temperature and pressure,
|
||||
* return the binary diffusion coefficent for components
|
||||
@@ -496,8 +425,8 @@ public:
|
||||
N2::gasEnthalpy(temperature, pN2);
|
||||
}
|
||||
else {
|
||||
Scalar pWater = fluidState.partialPressure(H2OIdx);
|
||||
Scalar pN2 = fluidState.partialPressure(N2Idx);
|
||||
Scalar pWater = fluidState.fugacity(H2OIdx);
|
||||
Scalar pN2 = fluidState.fugacity(N2Idx);
|
||||
|
||||
Scalar result = 0;
|
||||
result +=
|
||||
@@ -533,7 +462,7 @@ private:
|
||||
// \todo: proper citation
|
||||
Scalar rholH2O = H2O::liquidDensity(T, pl);
|
||||
Scalar clH2O = rholH2O/H2O::molarMass();
|
||||
|
||||
|
||||
// this assumes each nitrogen molecule displaces exactly one
|
||||
// water molecule in the liquid
|
||||
return
|
||||
@@ -542,69 +471,10 @@ private:
|
||||
xlN2*N2::molarMass());
|
||||
}
|
||||
|
||||
// defect of gas density
|
||||
static inline Scalar gasDefect_(Scalar pH2O, Scalar T, Scalar pg, Scalar xgH2O, Scalar xgN2)
|
||||
{
|
||||
// this assumes that sum of the fugacities of the individual
|
||||
// components add up to the gas pressure
|
||||
return
|
||||
H2O::gasDensity(T, pH2O)/H2O::molarMass()*(xgH2O - 1) +
|
||||
xgH2O*N2::gasDensity(T, pg - pH2O)/N2::molarMass();
|
||||
}
|
||||
|
||||
static Scalar gasPhaseDensity_(Scalar T, Scalar pg, Scalar xgH2O, Scalar xgN2)
|
||||
{
|
||||
// assume ideal gas for the initial condition
|
||||
Scalar pH2O = pg*xgH2O;
|
||||
Scalar delta = 1e100;
|
||||
|
||||
// newton method. makes sure that the total pressure of the
|
||||
// gas phase is the sum of the individual gas phase fugacities
|
||||
// of the individual components
|
||||
for (int i = 0; i < 10; ++i) {
|
||||
Scalar f;
|
||||
Scalar df_dpH2O;
|
||||
|
||||
f = gasDefect_(pH2O, T, pg, xgH2O, xgN2);
|
||||
Scalar eps = (std::abs(pg) + 1)*1e-9;
|
||||
df_dpH2O = gasDefect_(pH2O + eps, T, pg, xgH2O, xgN2);
|
||||
df_dpH2O -= gasDefect_(pH2O - eps, T, pg, xgH2O, xgN2);
|
||||
df_dpH2O /= 2*eps;
|
||||
|
||||
delta = - f/df_dpH2O;
|
||||
pH2O += delta;
|
||||
|
||||
if (std::abs(delta) < 1e-10*pg) {
|
||||
return H2O::gasDensity(T, pH2O) + N2::gasDensity(T, pg - pH2O);
|
||||
}
|
||||
};
|
||||
DUNE_THROW(NumericalProblem,
|
||||
"Can not calculate the gas density at "
|
||||
<< " T=" << T
|
||||
<< " pg=" << pg
|
||||
<< " xgH2O=" << xgH2O
|
||||
<< " xgN2=" << xgN2
|
||||
<< "\n");
|
||||
return H2O::gasDensity(T, pg*xgH2O) + N2::gasDensity(T, pg*xgN2);
|
||||
};
|
||||
|
||||
template <class FluidState>
|
||||
static void checkConsistentGasDensity_(Scalar rhoToTest,
|
||||
Scalar pg,
|
||||
const FluidState &fluidState)
|
||||
{
|
||||
#ifndef NDEBUG
|
||||
Scalar rho = gasPhaseDensity_(fluidState.temperature(),
|
||||
pg,
|
||||
fluidState.moleFrac(gPhaseIdx, H2OIdx),
|
||||
fluidState.moleFrac(gPhaseIdx, N2Idx));
|
||||
static const Scalar eps = 1e-8;
|
||||
if (std::abs(rhoToTest/rho - 1.0) > eps) {
|
||||
DUNE_THROW(Dune::InvalidStateException,
|
||||
"Density of gas phase is inconsistent: rho1/rho2 = "
|
||||
<< rhoToTest/rho);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
};
|
||||
|
||||
} // end namepace
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: idealgas.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2009 by Andreas Lauser
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: settablephase.hh 3777 2010-06-24 06:46:46Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2010 by Andreas Lauser
|
||||
* Institute of Hydraulic Engineering *
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
// $Id: boxspatialparameters.hh 3783 2010-06-24 11:33:53Z bernd $
|
||||
// $Id$
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2010 by Andreas Lauser *
|
||||
* Institute of Hydraulic Engineering *
|
||||
@@ -39,10 +39,14 @@ class BoxSpatialParameters
|
||||
{
|
||||
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) Implementation;
|
||||
|
||||
enum {
|
||||
dimWorld = GridView::dimensionworld
|
||||
};
|
||||
|
||||
typedef typename GridView::template Codim<0>::Entity Element;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
|
||||
|
||||
typedef typename GridView::ctype CoordScalar;
|
||||
typedef Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld> Tensor;
|
||||
@@ -54,6 +58,41 @@ public:
|
||||
~BoxSpatialParameters()
|
||||
{}
|
||||
|
||||
/*!
|
||||
* \brief Returns the factor by which the volume of a sub control
|
||||
* volume needs to be multiplied in order to get cubic
|
||||
* meters.
|
||||
*
|
||||
* By default that's just 1.0
|
||||
*/
|
||||
Scalar extrusionFactorScv(const Element &element,
|
||||
const FVElementGeometry &fvElemGeom,
|
||||
int scvIdx) const
|
||||
{ return 1.0; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the factor by which the area of a sub control
|
||||
* volume face needs to be multiplied in order to get
|
||||
* square meters.
|
||||
*
|
||||
* By default it is the arithmetic mean of the extrusion factor of
|
||||
* the face's two sub-control volumes.
|
||||
*/
|
||||
Scalar extrusionFactorScvf(const Element &element,
|
||||
const FVElementGeometry &fvElemGeom,
|
||||
int scvfIdx) const
|
||||
{
|
||||
return
|
||||
0.5 *
|
||||
(asImp_().extrusionFactorScv(element,
|
||||
fvElemGeom,
|
||||
fvElemGeom.subContVolFace[scvfIdx].i)
|
||||
+
|
||||
asImp_().extrusionFactorScv(element,
|
||||
fvElemGeom,
|
||||
fvElemGeom.subContVolFace[scvfIdx].j));
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Averages the intrinsic permeability.
|
||||
*/
|
||||
@@ -82,6 +121,13 @@ public:
|
||||
for (int j = 0; j < dimWorld; ++j)
|
||||
result[i][j] = harmonicMean(K1[i][j], K2[i][j]);
|
||||
}
|
||||
|
||||
protected:
|
||||
Implementation &asImp_()
|
||||
{ return *static_cast<Implementation*>(this); }
|
||||
|
||||
const Implementation &asImp_() const
|
||||
{ return *static_cast<const Implementation*>(this); }
|
||||
};
|
||||
|
||||
} // namespace Dumux
|
||||
|
||||
Reference in New Issue
Block a user