diff --git a/dumux/material/binarycoefficients/fullermethod.hh b/dumux/material/binarycoefficients/fullermethod.hh index b7895d135..decea899c 100644 --- a/dumux/material/binarycoefficients/fullermethod.hh +++ b/dumux/material/binarycoefficients/fullermethod.hh @@ -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 * diff --git a/dumux/material/binarycoefficients/h2o_n2.hh b/dumux/material/binarycoefficients/h2o_n2.hh index aa152aee9..682abd2a9 100644 --- a/dumux/material/binarycoefficients/h2o_n2.hh +++ b/dumux/material/binarycoefficients/h2o_n2.hh @@ -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 * diff --git a/dumux/material/binarycoefficients/henryiapws.hh b/dumux/material/binarycoefficients/henryiapws.hh index f0baa26e8..df4c30afd 100644 --- a/dumux/material/binarycoefficients/henryiapws.hh +++ b/dumux/material/binarycoefficients/henryiapws.hh @@ -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 * diff --git a/dumux/material/components/component.hh b/dumux/material/components/component.hh index 466d3737e..9ae42e9ed 100644 --- a/dumux/material/components/component.hh +++ b/dumux/material/components/component.hh @@ -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 - namespace Dumux { @@ -33,6 +31,10 @@ template 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. */ diff --git a/dumux/material/components/h2.hh b/dumux/material/components/h2.hh index c736b782e..31d707bb7 100644 --- a/dumux/material/components/h2.hh +++ b/dumux/material/components/h2.hh @@ -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); diff --git a/dumux/material/components/h2o.hh b/dumux/material/components/h2o.hh index a197c7285..53dd670ad 100644 --- a/dumux/material/components/h2o.hh +++ b/dumux/material/components/h2o.hh @@ -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 * diff --git a/dumux/material/components/iapws/common.hh b/dumux/material/components/iapws/common.hh index 659b7a59b..991cca533 100644 --- a/dumux/material/components/iapws/common.hh +++ b/dumux/material/components/iapws/common.hh @@ -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 * diff --git a/dumux/material/components/iapws/region1.hh b/dumux/material/components/iapws/region1.hh index 528db976f..cbcb8be2f 100644 --- a/dumux/material/components/iapws/region1.hh +++ b/dumux/material/components/iapws/region1.hh @@ -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 * diff --git a/dumux/material/components/iapws/region2.hh b/dumux/material/components/iapws/region2.hh index c8fb50d44..f04599881 100644 --- a/dumux/material/components/iapws/region2.hh +++ b/dumux/material/components/iapws/region2.hh @@ -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 * diff --git a/dumux/material/components/iapws/region4.hh b/dumux/material/components/iapws/region4.hh index e8c5066c5..c527846b3 100644 --- a/dumux/material/components/iapws/region4.hh +++ b/dumux/material/components/iapws/region4.hh @@ -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 * diff --git a/dumux/material/components/n2.hh b/dumux/material/components/n2.hh index 1b9ad9662..c40943810 100644 --- a/dumux/material/components/n2.hh +++ b/dumux/material/components/n2.hh @@ -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 * diff --git a/dumux/material/components/o2.hh b/dumux/material/components/o2.hh index 7fa102bf5..6cef5df49 100644 --- a/dumux/material/components/o2.hh +++ b/dumux/material/components/o2.hh @@ -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); diff --git a/dumux/material/components/simplednapl.hh b/dumux/material/components/simplednapl.hh index a00a7d34f..2ec18d3b3 100644 --- a/dumux/material/components/simplednapl.hh +++ b/dumux/material/components/simplednapl.hh @@ -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 * diff --git a/dumux/material/components/simpleh2o.hh b/dumux/material/components/simpleh2o.hh index ad8fd2aa6..601b20993 100644 --- a/dumux/material/components/simpleh2o.hh +++ b/dumux/material/components/simpleh2o.hh @@ -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 * diff --git a/dumux/material/components/tabulatedcomponent.hh b/dumux/material/components/tabulatedcomponent.hh index 588cf6d14..7cf871f74 100644 --- a/dumux/material/components/tabulatedcomponent.hh +++ b/dumux/material/components/tabulatedcomponent.hh @@ -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 * diff --git a/dumux/material/components/unit.hh b/dumux/material/components/unit.hh index 31b3d4afb..333b290fd 100644 --- a/dumux/material/components/unit.hh +++ b/dumux/material/components/unit.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh index 5cd1d5d53..2bf6fa2f0 100644 --- a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh +++ b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh b/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh index 3bf6a5839..ae4bf6cdd 100644 --- a/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh +++ b/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh b/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh index 96f4c92eb..c8369dbed 100644 --- a/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh +++ b/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh b/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh index a1a7049e6..4ef065595 100644 --- a/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh +++ b/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh index f8b6dd43a..7e96ef2b3 100644 --- a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh +++ b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/linearmaterialparams.hh b/dumux/material/fluidmatrixinteractions/2p/linearmaterialparams.hh index f3a9fe506..b0c95985a 100644 --- a/dumux/material/fluidmatrixinteractions/2p/linearmaterialparams.hh +++ b/dumux/material/fluidmatrixinteractions/2p/linearmaterialparams.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh index 51d7952be..e5b4b664b 100644 --- a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh +++ b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh index 767a47fb0..c65b2869e 100644 --- a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh +++ b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh index 072df298d..80233098b 100644 --- a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh +++ b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh index bbac9ad8d..c65ff9f7a 100644 --- a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh +++ b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh index 45102ffe6..88603f05c 100644 --- a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh +++ b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh @@ -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 * diff --git a/dumux/material/fluidmatrixinteractions/2p/vangenuchtenparams.hh b/dumux/material/fluidmatrixinteractions/2p/vangenuchtenparams.hh index 25181813b..0812da791 100644 --- a/dumux/material/fluidmatrixinteractions/2p/vangenuchtenparams.hh +++ b/dumux/material/fluidmatrixinteractions/2p/vangenuchtenparams.hh @@ -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 * diff --git a/dumux/material/fluidstate.hh b/dumux/material/fluidstate.hh index e496438f9..1d9dc2599 100644 --- a/dumux/material/fluidstate.hh +++ b/dumux/material/fluidstate.hh @@ -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. diff --git a/dumux/material/fluidsystems/2p_system.hh b/dumux/material/fluidsystems/2p_system.hh index e1dab713b..eea985f1e 100644 --- a/dumux/material/fluidsystems/2p_system.hh +++ b/dumux/material/fluidsystems/2p_system.hh @@ -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 * diff --git a/dumux/material/fluidsystems/h2o_n2_system.hh b/dumux/material/fluidsystems/h2o_n2_system.hh index a4154846c..f1eb18b80 100644 --- a/dumux/material/fluidsystems/h2o_n2_system.hh +++ b/dumux/material/fluidsystems/h2o_n2_system.hh @@ -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 - 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 - 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 - 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 diff --git a/dumux/material/idealgas.hh b/dumux/material/idealgas.hh index 9bc012277..5592c17b5 100644 --- a/dumux/material/idealgas.hh +++ b/dumux/material/idealgas.hh @@ -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 * diff --git a/dumux/material/settablephase.hh b/dumux/material/settablephase.hh index 29320b578..621ce608c 100644 --- a/dumux/material/settablephase.hh +++ b/dumux/material/settablephase.hh @@ -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 * diff --git a/dumux/material/spatialparameters/boxspatialparameters.hh b/dumux/material/spatialparameters/boxspatialparameters.hh index f15cdd6b3..3a1a4b97b 100644 --- a/dumux/material/spatialparameters/boxspatialparameters.hh +++ b/dumux/material/spatialparameters/boxspatialparameters.hh @@ -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 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(this); } + + const Implementation &asImp_() const + { return *static_cast(this); } }; } // namespace Dumux