// -*- 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
* \ingroup Fluidsystems
* \brief A fluid system with water and gas as phases and brine and CO2
* as components.
*/
#ifndef DUMUX_BIOMIN_FLUID_SYSTEM_HH
#define DUMUX_BIOMIN_FLUID_SYSTEM_HH
#include <dumux/common/parameters.hh>
#include <dumux/material/idealgas.hh>
#include <dumux/material/fluidsystems/base.hh>
#include <dumux/material/components/brine.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/co2.hh>
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/calciumion.hh>
#include "../components/urea.hh"
#include <dumux/material/binarycoefficients/brine_co2.hh>
namespace Dumux {
namespace FluidSystems {
/*!
* \ingroup Fluidsystems
* \brief A compositional fluid with brine and carbon dioxide as
* components in both, the liquid and the gas (supercritical) phase,
* additional biomineralisation components (Ca and Urea) in the liquid phase
*
* This class provides acess to the Bio fluid system when no property system is used.
* For Dumux users, using BioMinFluid<TypeTag> and the documentation therein is
* recommended.
*
* The user can provide their own material table for co2 properties.
* This fluidsystem is initialized as default with the tabulated version of
* water of the IAPWS-formulation, and the tabularized adapter to transfer
* this into brine.
* In the non-TypeTagged version, salinity information has to be provided with
* the init() methods.
*/
template <class Scalar,
class CO2Table,
class H2OType = Dumux::Components::TabulatedComponent<Dumux::Components::H2O<Scalar>> >
class BioMin
: public Base<Scalar, BioMin<Scalar, CO2Table, H2OType> >
{
using ThisType = BioMin<Scalar, H2OType>;
using Base = Dumux::FluidSystems::Base<Scalar, ThisType>;
using Brine_CO2 = BinaryCoeff::Brine_CO2<Scalar, CO2Table>;
using IdealGas = Dumux::IdealGas<Scalar>;
public:
using CO2 = Components::CO2<Scalar, CO2Table>;
using H2O = H2OType;
using Ca = Components::CalciumIon<Scalar>;
using Urea = Components::Urea<Scalar>;
using Brine = Components::Brine<Scalar, H2O>;
// the type of parameter cache objects. this fluid system does not
// cache anything, so it uses Dumux::NullParameterCache
using ParameterCache = Dumux::NullParameterCache;
/****************************************
* Fluid phase related static parameters
****************************************/
static constexpr int numPhases = 2; // liquid and gas phases
static constexpr int liquidPhaseIdx = 0; // index of the liquid phase
static constexpr int gasPhaseIdx = 1; // index of the gas phase
static constexpr int phase0Idx = liquidPhaseIdx; // index of the first phase
static constexpr int phase1Idx = gasPhaseIdx; // index of the second phase
/*!
* \brief Return the human readable name of a fluid phase
*
* \param phaseIdx The index of the fluid phase to consider
*/
static std::string phaseName(int phaseIdx)
{
static std::string name[] =
{
"liquid",
"gas"
};
assert(0 <= phaseIdx && phaseIdx < numPhases);
return name[phaseIdx];
}
/*!
* \brief Returns whether the fluids are miscible
*/
static constexpr bool isMiscible()
{ return true; }
/*!
* \brief Return whether a phase is liquid
*
* \param phaseIdx The index of the fluid phase to consider
*/
static constexpr bool isLiquid(int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
return phaseIdx != gasPhaseIdx;
}
/*!
* \brief Returns true if and only if a fluid phase is assumed to
* be an ideal mixture.
*
* We define an ideal mixture as a fluid phase where the fugacity
* coefficients of all components times the pressure of the phase
* are independent on the fluid composition. This assumption is true
* if Henry's law and Raoult's law apply. If you are unsure what
* this function should return, it is safe to return false. The
* only damage done will be (slightly) increased computation times
* in some cases.
*
* \param phaseIdx The index of the fluid phase to consider
*/
static constexpr bool isIdealMixture(int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
return true;
}
/*!
* \brief Returns true if and only if a fluid phase is assumed to
* be compressible.
*
* Compressible means that the partial derivative of the density
* to the fluid pressure is always larger than zero.
*
* \param phaseIdx The index of the fluid phase to consider
*/
static constexpr bool isCompressible(int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
return true;
}
/*!
* \brief Returns true if and only if a fluid phase is assumed to
* be an ideal gas.
*
* \param phaseIdx The index of the fluid phase to consider
*/
static bool isIdealGas(int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
// let the fluids decide
if (phaseIdx == gasPhaseIdx)
return H2O::gasIsIdeal() && CO2::gasIsIdeal();
return false; // not a gas
}
/****************************************
* Component related static parameters
****************************************/
static constexpr int numComponents = 4; // H2O/brine, CO2, Ca, urea
static constexpr int H2OIdx = 0;
static constexpr int CO2Idx = 1;
static constexpr int CaIdx = 2;
static constexpr int UreaIdx = 3;
static constexpr int BrineIdx = H2OIdx;
static constexpr int comp0Idx = BrineIdx;
static constexpr int comp1Idx = CO2Idx;
/*!
* \brief Return the human readable name of a component
*
* \param compIdx The index of the component to consider
*/
static std::string componentName(int compIdx)
{
switch (compIdx) {
case BrineIdx: return Brine::name();
case CO2Idx: return "CO2";
case CaIdx: return Ca::name();
case UreaIdx: return Urea::name();
default: DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
};
}
/*!
* \brief Return the molar mass of a component in \f$\mathrm{[kg/mol]}\f$.
*
* \param compIdx The index of the component to consider
*/
static Scalar molarMass(int compIdx)
{
switch (compIdx) {
case H2OIdx: return H2O::molarMass();
// actually, the molar mass of brine is only needed for diffusion
// but since solutes are accounted for seperately
// only the molar mass of water is returned.
case CO2Idx: return CO2::molarMass();
case CaIdx: return Ca::molarMass();
case UreaIdx: return Urea::molarMass();
default: DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
};
}
/****************************************
* thermodynamic relations
****************************************/
static void init()
{
init(/*startTemp=*/295.15, /*endTemp=*/305.15, /*tempSteps=*/10,
/*startPressure=*/1e4, /*endPressure=*/1e6, /*pressureSteps=*/2000);
}
static void init(Scalar startTemp, Scalar endTemp, int tempSteps,
Scalar startPressure, Scalar endPressure, int pressureSteps)
{
std::cout << "Initializing tables for the pure-water properties.\n";
H2O::init(startTemp, endTemp, tempSteps,
startPressure, endPressure, pressureSteps);
}
/*!
* \brief Given a phase's composition, temperature, pressure, and
* the partial pressures of all components, return its
* density \f$\mathrm{[kg/m^3]}\f$.
*
* \param fluidState The fluid state
* \param phaseIdx The index of the phase
*/
template <class FluidState>
static Scalar density(const FluidState &fluidState,
const ParameterCache ¶mCache,
int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
switch (phaseIdx) {
// assume pure brine for the liquid phase.
case liquidPhaseIdx:
return Brine::liquidDensity(temperature, pressure); //TODO major assumption in favor of runtime! Not considering density effect of dissolved calcium
// assume pure CO2 for the gas phase.
case gasPhaseIdx:
return CO2::gasDensity(temperature, pressure);
default:
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); break;
}
}
/*!
* \brief The molar density \f$\rho_{mol,\alpha}\f$
* of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
*
* The molar density is defined by the
* mass density \f$\rho_\alpha\f$ and the mean molar mass \f$\overline M_\alpha\f$:
*
* \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
*/
template <class FluidState>
static Scalar molarDensity(const FluidState &fluidState,
const ParameterCache ¶mCache,
int phaseIdx)
{
const Scalar temperature = fluidState.temperature(phaseIdx);
const Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == liquidPhaseIdx)
{
return density(fluidState, paramCache, phaseIdx)
/ Brine::molarMass();
}
else if (phaseIdx == gasPhaseIdx)
{
// for the gas phase assume an ideal gas
return CO2::gasMolarDensity(temperature, pressure);
}
else
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
/*!
* \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*
* Equation given in: - Batzle & Wang (1992)
* - cited by: Bachu & Adams (2002)
* "Equations of State for basin geofluids"
*/
template <class FluidState>
static Scalar viscosity(const FluidState &fluidState,
const ParameterCache ¶mCache,
int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
// assume brine with viscosity effect of Ca for the liquid phase.
if (phaseIdx == liquidPhaseIdx)
return Brine::liquidViscosity(temperature, pressure);
// assume pure CO2 for the gas phase.
else if (phaseIdx == gasPhaseIdx)
return CO2::gasViscosity(temperature, pressure);
else
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
/*!
* \brief Returns the fugacity coefficient [Pa] of a component in a
* phase.
*
* The fugacity coefficient \f$\phi^\kappa_\alpha\f$ of a
* component \f$\kappa\f$ for a fluid phase \f$\alpha\f$ defines
* the fugacity \f$f^\kappa_\alpha\f$ by the equation
*
* \f[
f^\kappa_\alpha := \phi^\kappa_\alpha x^\kappa_\alpha p_\alpha\;.
\f]
*
* The fugacity itself is just an other way to express the
* chemical potential \f$\zeta^\kappa_\alpha\f$ of the component:
*
* \f[
f^\kappa_\alpha := \exp\left\{\frac{\zeta^\kappa_\alpha}{k_B T_\alpha} \right\}
\f]
* where \f$k_B\f$ is Boltzmann's constant.
*/
template <class FluidState>
static Scalar fugacityCoefficient(const FluidState &fluidState,
const ParameterCache ¶mCache,
int phaseIdx,
int compIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
if (phaseIdx == gasPhaseIdx)
// use the fugacity coefficients of an ideal gas. the
// actual value of the fugacity is not relevant, as long
// as the relative fluid compositions are observed,
return 1.0;
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
Scalar salinity = Brine::salinity(); // 0.1; //TODO major assumption in favor of runtime!
//function is actually designed for use with NaCl not Ca.
//Theoretically it should be: fluidState.massFraction(liquidPhaseIdx, CaIdx);
assert(temperature > 0);
assert(pressure > 0);
// calulate the equilibrium composition for the given
// temperature and pressure.
Scalar xwH2O, xnH2O;
Scalar xwCO2, xnCO2;
Brine_CO2::calculateMoleFractions(temperature,
pressure,
salinity,
/*knowgasPhaseIdx=*/-1,
xwCO2,
xnH2O);
// normalize the phase compositions
using std::min;
using std::max;
xwCO2 = max(0.0, min(1.0, xwCO2));
xnH2O = max(0.0, min(1.0, xnH2O));
xwH2O = 1.0 - xwCO2;
xnCO2 = 1.0 - xnH2O;
if (compIdx == BrineIdx)
{
Scalar phigH2O = 1.0;
return phigH2O * xnH2O / xwH2O;
}
if (compIdx == CO2Idx)
{
Scalar phigCO2 = 1.0;
return phigCO2 * xnCO2 / xwCO2;
}
else
return 1/pressure; //all other components stay in the liquid phase
}
/*!
* \brief Given the phase compositions, return the binary
* diffusion coefficient \f$\mathrm{[m^2/s]}\f$ of two components in a phase.
* \param fluidState An arbitrary fluid state
* \param phaseIdx The index of the fluid phase to consider
* \param compIIdx Index of the component i
* \param compJIdx Index of the component j
*/
template <class FluidState>
static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
const ParameterCache ¶mCache,
int phaseIdx,
int compIIdx,
int compJIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIIdx && compIIdx < numComponents);
assert(0 <= compJIdx && compJIdx < numComponents);
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == liquidPhaseIdx)
{
assert(compIIdx == H2OIdx);
if (compJIdx == CO2Idx)
return Brine_CO2::liquidDiffCoeff(temperature, pressure);
else if (compJIdx < numComponents) //Calcium and urea
return 1.587e-9; //[m²/s] educated guess, value for NaCl from J. Phys. D: Appl. Phys. 40 (2007) 2769-2776
else
DUNE_THROW(Dune::NotImplemented, "Binary difussion coefficient : Incorrect compIdx");
}
else
{
assert(phaseIdx == gasPhaseIdx);
assert(compIIdx == CO2Idx);
if (compJIdx == H2OIdx)
return Brine_CO2::gasDiffCoeff(temperature, pressure);
else if (compJIdx <numComponents) //Calcium and urea will stay in brine, no gaseous calcium or urea!
return 0.0;
else
DUNE_THROW(Dune::NotImplemented, "Binary difussion coefficient : Incorrect compIdx");
}
};
/*!
* \brief Given a phase's composition, temperature and pressure,
* return its specific enthalpy \f$\mathrm{[J/kg]}\f$.
* \param fluidState An arbitrary fluid state
* \param phaseIdx The index of the fluid phase to consider
*
* See:
* Class 2001:
* Theorie und numerische Modellierung nichtisothermer Mehrphasenprozesse in NAPL-kontaminierten porösen Medien
* Chapter 2.1.13 Innere Energie, Wäremekapazität, Enthalpie \cite A3:class:2001 <BR>
*
* Formula (2.42):
* the specific enthalpy of a gasphase result from the sum of (enthalpies*mass fraction) of the components
*
*/
template <class FluidState>
static Scalar enthalpy(const FluidState &fluidState,
const ParameterCache ¶mCache,
int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
Scalar temperature = fluidState.temperature(phaseIdx);
Scalar pressure = fluidState.pressure(phaseIdx);
if (phaseIdx == liquidPhaseIdx)
{
// assume pure brine for the liquid phase.
return Brine::liquidEnthalpy(temperature, pressure);
}
else
{
// assume pure CO2 for the gas phase.
return CO2::gasEnthalpy(temperature, pressure);
}
};
};
} // end namespace FluidSystems
} // end namespace Dumux
#endif