Added lots of stuff originally from the samcode repository.

This commit is contained in:
Atgeirr Flø Rasmussen 2010-11-12 13:18:27 +01:00
commit 45bf24c9b1
16 changed files with 1634 additions and 0 deletions

View File

@ -0,0 +1,165 @@
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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 3 of the License, or
(at your option) any later version.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILCOMPONENT_HEADER_INCLUDED
#define OPM_BLACKOILCOMPONENT_HEADER_INCLUDED
#include <dune/common/stdstreams.hh>
namespace Dumux
{
/*!
* \ingroup Components
* \brief
* A component class for the black oil model, intended to be used
* for all three components.
*
* \tparam Scalar The type used for scalar values
*/
template <class Scalar>
class BlackoilComponent
{
public:
/*!
* \brief A human readable name for the component.
*/
static const char *name()
{
return "BlackoilComponent";
}
/*!
* \brief The molar mass in [kg] of the component.
*/
static Scalar molarMass()
{ DUNE_THROW(Dune::NotImplemented, "Component::molarMass()"); }
/*!
* \brief Returns the critical temperature in [K] of the component.
*/
static Scalar criticalTemperature()
{ DUNE_THROW(Dune::NotImplemented, "Component::criticalTemperature()"); }
/*!
* \brief Returns the critical pressure in [Pa] of the component.
*/
static Scalar criticalPressure()
{ DUNE_THROW(Dune::NotImplemented, "Component::criticalPressure()"); }
/*!
* \brief Returns the temperature in [K] at the component's triple point.
*/
static Scalar tripleTemperature()
{ DUNE_THROW(Dune::NotImplemented, "Component::tripleTemperature()"); }
/*!
* \brief Returns the pressure in [Pa] at the component's triple point.
*/
static Scalar triplePressure()
{ DUNE_THROW(Dune::NotImplemented, "Component::triplePressure()"); }
/*!
* \brief The vapor pressure in [Pa] of the component at a given
* temperature in [K].
*
* \param T temperature of the component in [K]
*/
static Scalar vaporPressure(Scalar T)
{ DUNE_THROW(Dune::NotImplemented, "Component::vaporPressure()"); }
/*!
* \brief The density in [kg/m^3] of the component at a given pressure in [Pa] and temperature in [K].
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::density()"); }
/*!
* \brief The density [kg/m^3] of the liquid component at a given pressure in [Pa] and temperature in [K].
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::density()"); }
/*!
* \brief Specific enthalpy [J/kg] of the pure component in gas.
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::gasEnthalpy()"); }
/*!
* \brief Specific enthalpy [J/kg] of the pure component in liquid.
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::liquidEnthalpy()"); }
/*!
* \brief Specific internal energy [J/kg] of the pure component in gas.
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::gasInternalEnergy()"); }
/*!
* \brief Specific internal energy [J/kg] of pure the pure component in liquid.
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::liquidInternalEnergy()"); }
/*!
* \brief The dynamic viscosity [Pa*s] of the pure component at a given pressure in [Pa] and temperature in [K].
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::gasViscosity()"); }
/*!
* \brief The dynamic liquid viscosity [Pa*s] of the pure component.
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::liquidViscosity()"); }
};
} // end namepace
#endif // OPM_BLACKOILCOMPONENT_HEADER_INCLUDED

View File

@ -0,0 +1,19 @@
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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 3 of the License, or
(at your option) any later version.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/

View File

@ -0,0 +1,66 @@
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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 3 of the License, or
(at your option) any later version.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILPVT_HEADER_INCLUDED
#define OPM_BLACKOILPVT_HEADER_INCLUDED
#include <string>
#include <boost/scoped_ptr.hpp>
#include "MiscibilityProps.hpp"
namespace Opm
{
class BlackoilPVT
{
public:
typedef MiscibilityProps::surfvol_t surfvol_t;
enum PhaseIndex { Aqua = 0, Vapour = 1, Liquid = 2 };
void init(const Dune::EclipseGridParser& ep);
double getViscosity(double press, const surfvol_t& surfvol,
PhaseNames phase) const;
surfvol_t getMobilities(double press, const surfvol_t& sat, const surfvol_t& surfvol) const;
surfvol_t surfaceDensities() const;
double B (double press, const surfvol_t& surfvol,
PhaseNames phase) const;
double dBdp(double press, const surfvol_t& surfvol,
PhaseNames phase) const;
double R (double press, const surfvol_t& surfvol,
PhaseNames phase) const;
double dRdp(double press, const surfvol_t& surfvol,
PhaseNames phase) const;
private:
int region_number_;
const MiscibilityProps& propsForPhase(PhaseNames phase) const;
boost::scoped_ptr<MiscibilityProps> water_props_;
boost::scoped_ptr<MiscibilityProps> oil_props_;
boost::scoped_ptr<MiscibilityProps> gas_props_;
surfvol_t densities_;
};
}
#endif // OPM_BLACKOILPVT_HEADER_INCLUDED

View File

@ -0,0 +1,242 @@
//===========================================================================
//
// File: FluidMiscibilityThreePhase.cpp
//
// Created: Wed Feb 10 09:25:57 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
#include "FluidMiscibilityThreePhase.hpp"
#include <dune/common/EclipseGridParser.hpp>
#include "MiscibilityDead.hpp"
#include "MiscibilityLiveOil.hpp"
#include "MiscibilityLiveGas.hpp"
#include "MiscibilityWater.hpp"
#include <dune/common/ErrorMacros.hpp>
#include <dune/common/linInt.hpp>
using namespace Dune;
namespace samcode
{
void FluidMiscibilityThreePhase::init(const std::string& pvt_filename)
{
typedef std::vector<std::vector<std::vector<double> > > table_t;
region_number_ = 0;
Dune::EclipseGridParser eclipse_props(pvt_filename);
// Surface densities. Accounting for different orders in eclipse and our code.
if (eclipse_props.hasField("DENSITY")) {
const int region_number = 0;
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
const std::vector<std::vector<double> > d_tmp =
eclipse_props.getDENSITY().densities_;
densities_[Aqua] = d_tmp[region_number][ECL_water];
densities_[Liquid] = d_tmp[region_number][ECL_oil];
densities_[Vapour] = d_tmp[region_number][ECL_gas];
} else {
THROW("DENSITY not defined");
}
// Water PVT
if (eclipse_props.hasField("PVTW")) {
water_props_.reset(new MiscibilityWater(eclipse_props.getPVTW().pvtw_));
} else {
water_props_.reset(new MiscibilityWater(3e-4)); // Default is 0.3 cP.
}
// Oil PVT
if (eclipse_props.hasField("PVDO")) {
oil_props_.reset(new MiscibilityDead(eclipse_props.getPVDO().pvdo_));
} else if (eclipse_props.hasField("PVTO")) {
oil_props_.reset(new MiscibilityLiveOil(eclipse_props.getPVTO().pvto_));
} else {
THROW("File " << pvt_filename << " is missing PVDO and PVTO\n");
}
// Gas PVT
if (eclipse_props.hasField("PVDG")) {
gas_props_.reset(new MiscibilityDead(eclipse_props.getPVDG().pvdg_));
} else if (eclipse_props.hasField("PVTG")) {
gas_props_.reset(new MiscibilityLiveGas(eclipse_props.getPVTG().pvtg_));
} else {
THROW("File " << pvt_filename << " is missing PVDG and PVTG\n");
}
//SWOF
if (eclipse_props.hasField("SWOF")) {
swof_ = eclipse_props.getSWOF().swof_;
}
//SGOF
if (eclipse_props.hasField("SGOF")) {
sgof_ = eclipse_props.getSGOF().sgof_;
}
}
FluidMiscibilityThreePhase::surfvol_t FluidMiscibilityThreePhase::surfaceDensities() const
{
return densities_;
}
double FluidMiscibilityThreePhase::getViscosity(double press, const surfvol_t& surfvol, PhaseNames phase) const
{
return propsForPhase(phase).getViscosity(region_number_, press, surfvol);
}
FluidMiscibilityThreePhase::surfvol_t
FluidMiscibilityThreePhase::getMobilities(double press, const surfvol_t& sat,
const surfvol_t& surfvol) const
{
if (swof_.empty() || sgof_.empty()) {
THROW("The SWOF and SGOF keywords were not given, cannot compute mobilities. Try tracer_flow=true.");
}
surfvol_t mobilities;
double sw = sat[Aqua];
double sg = sat[Vapour];
mobilities[Aqua] = krw(sw)/getViscosity(press, surfvol, Aqua);
mobilities[Liquid] = krow(sw)*krog(sg)/getViscosity(press, surfvol, Liquid);
mobilities[Vapour] = krg(sg)/getViscosity(press, surfvol, Vapour);
return mobilities;
}
double FluidMiscibilityThreePhase::B(double press, const surfvol_t& surfvol, PhaseNames phase) const
{
return propsForPhase(phase).B(region_number_, press, surfvol);
}
double FluidMiscibilityThreePhase::dBdp(double press, const surfvol_t& surfvol, PhaseNames phase) const
{
return propsForPhase(phase).dBdp(region_number_, press, surfvol);
}
double FluidMiscibilityThreePhase::R(double press, const surfvol_t& surfvol, PhaseNames phase) const
{
return propsForPhase(phase).R(region_number_, press, surfvol);
}
double FluidMiscibilityThreePhase::dRdp(double press, const surfvol_t& surfvol, PhaseNames phase) const
{
return propsForPhase(phase).dRdp(region_number_, press, surfvol);
}
const MiscibilityProps& FluidMiscibilityThreePhase::propsForPhase(PhaseNames phase) const
{
switch (phase) {
case Aqua:
return *water_props_;
case Liquid:
return *oil_props_;
case Vapour:
return *gas_props_;
default:
THROW("Unknown phase accessed: " << phase);
}
}
// Stone's first model(Modified)
double FluidMiscibilityThreePhase::kro(double sw, double sg) const
{
double so = 1.0-sw-sg;
double fw = krow(sw);
double fg = krog(sg);
return so*fw*fg;
}
// kro partial derivative Sw
double FluidMiscibilityThreePhase::dkro_dsw(double sw, double sg) const
{
double so = 1.0-sw-sg;
double fw = krow(sw);
double fg = krog(sg);
double fwder = dkrow_dsw(sw);
return fg*(-fw + so*fwder);
}
// kro partial derivative Sg
double FluidMiscibilityThreePhase::dkro_dsg(double sw, double sg) const
{
double so = 1.0-sw-sg;
double fw = krow(sw);
double fg = krog(sg);
double fgder = dkrog_dsg(sg);
return fw*(-fg + so*fgder);
}
// kro partial derivatives
void FluidMiscibilityThreePhase::dkro(double sw, double sg, double& dkro_dsw,
double& dkro_dsg) const
{
double so = 1.0-sw-sg;
double fw = krow(sw);
double fg = krog(sg);
double fwder = dkrow_dsw(sw);
double fgder = dkrog_dsg(sg);
dkro_dsw = fg*(-fw + so*fwder);
dkro_dsg = fw*(-fg + so*fgder);
}
// Water relative permeability
double FluidMiscibilityThreePhase::krw (double sw) const
{
return linearInterpolation(swof_[region_number_][0],
swof_[region_number_][1], sw);
}
// Oil relative permeability
double FluidMiscibilityThreePhase::krow (double sw) const
{
return linearInterpolation(swof_[region_number_][0],
swof_[region_number_][2], sw);
}
// krow derivative
double FluidMiscibilityThreePhase::dkrow_dsw(double sw) const
{
return linearInterpolDerivative(swof_[region_number_][0],
swof_[region_number_][2], sw);
}
// Water-oil capillary pressure
double FluidMiscibilityThreePhase::Pcow (double sw) const
{
return linearInterpolation(swof_[region_number_][0],
swof_[region_number_][3], sw);
}
// Gas relative permeability
double FluidMiscibilityThreePhase::krg (double sg) const
{
return linearInterpolation(sgof_[region_number_][0],
sgof_[region_number_][1], sg);
}
// Oil relative permeability
double FluidMiscibilityThreePhase::krog (double sg) const
{
return linearInterpolation(sgof_[region_number_][0],
sgof_[region_number_][2], sg);
}
// krog derivative
double FluidMiscibilityThreePhase::dkrog_dsg(double sg) const
{
return linearInterpolDerivative(sgof_[region_number_][0],
sgof_[region_number_][2], sg);
}
// Oil-gas capillary pressure
double FluidMiscibilityThreePhase::Pcog(double sg) const
{
return linearInterpolation(sgof_[region_number_][0],
sgof_[region_number_][3], sg);
}
} // namespace samcode

View File

@ -0,0 +1,78 @@
//===========================================================================
//
// File: FluidMiscibilityThreePhase.hpp
//
// Created: Thu Feb 11 10:35:05 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
#ifndef SINTEF_FLUIDMISCIBILITYTHREEPHASE_HEADER
#define SINTEF_FLUIDMISCIBILITYTHREEPHASE_HEADER
/** Temporary class for testing Miscibility* classes
* Detailed description.
*/
#include <string>
#include <boost/scoped_ptr.hpp>
#include "MiscibilityProps.hpp"
namespace samcode
{
class FluidMiscibilityThreePhase
{
public:
typedef MiscibilityProps::surfvol_t surfvol_t;
enum PhaseNames { Aqua = 0, Liquid = 1, Vapour = 2 };
FluidMiscibilityThreePhase(){}
~FluidMiscibilityThreePhase(){}
void init(const std::string& pvt_filename);
double getViscosity(double press, const surfvol_t& surfvol,
PhaseNames phase) const;
surfvol_t getMobilities(double press, const surfvol_t& sat, const surfvol_t& surfvol) const;
surfvol_t surfaceDensities() const;
double B (double press, const surfvol_t& surfvol,
PhaseNames phase) const;
double dBdp(double press, const surfvol_t& surfvol,
PhaseNames phase) const;
double R (double press, const surfvol_t& surfvol,
PhaseNames phase) const;
double dRdp(double press, const surfvol_t& surfvol,
PhaseNames phase) const;
private:
int region_number_;
const MiscibilityProps& propsForPhase(PhaseNames phase) const;
double kro(double sw, double sg) const; // Stone's first model(Modified)
double dkro_dsw(double sw, double sg) const; // kro partial derivative dSw
double dkro_dsg(double sw, double sg) const; // kro partial derivative dSg
void dkro (double sw, double sg, // kro partial derivatives. dSw and dSg
double& dkro_dsw, double& dkro_dsg) const;
double krw (double sw) const; // Water relative permeability
double krow (double sw) const; // Oil relative permeability
double dkrow_dsw (double sw) const; // krow derivative
double Pcow (double sw) const; // Water-oil capillary pressure
double krg (double sg) const; // Gas relative permeability
double krog (double sg) const; // Oil relative permeability
double dkrog_dsg(double sg) const; // krog derivative
double Pcog (double sg) const; // Oil-gas capillary pressure
boost::scoped_ptr<MiscibilityProps> water_props_;
boost::scoped_ptr<MiscibilityProps> oil_props_;
boost::scoped_ptr<MiscibilityProps> gas_props_;
surfvol_t densities_;
std::vector<std::vector<std::vector<double> > > sgof_; // Gas/Oil saturation
std::vector<std::vector<std::vector<double> > > swof_; // Water/Oil saturation
};
}
#endif // SINTEF_FLUIDMISCIBILITYTHREEPHASE_HEADER

View File

@ -0,0 +1,244 @@
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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 3 of the License, or
(at your option) any later version.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED
#define OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED
namespace Opm
{
/*!
* \brief Calcultes the phase state from the primary variables in the
* blackoil model.
*/
class FluidStateBlackoil
{
public:
enum { numPhases = 3 };
enum { numComponents = 3};
typedef int FluidSystem; // TODO
typedef int PrimaryVars; // TODO
typedef double Scalar;
typedef FluidMatrixInteractionBlackoil<Scalar> MaterialLaw;
typedef typename MaterialLaw::Params MaterialLawParams;
/*!
* \brief Update the phase state from the primary variables.
*/
void update(const PrimaryVariables &primaryVars,
const MaterialLawParams &pcParams,
Scalar temperature)
{
// calculate the temperature
temperature_ = temperature;
// calculate the saturations
saturation_[numPhases - 1] = 1.0;
for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
saturation_[phaseIdx] = primaryVars[S0Idx + phaseIdx];
saturation_[numPhases - 1] -= saturation_[phaseIdx];
}
// let the material law calculate the capillary pressures
MaterialLaw::pC(phasePressure_,
pcParams,
saturation_,
temperature_);
// convert to phase pressures
Scalar sumPc = 0;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
sumPc += phasePressure_[phaseIdx];
phasePressure_[phaseIdx] = primaryVars[p0Idx] - sumPc;
}
updateComposition_(primaryVars);
}
private:
void updateComposition_(const PrimaryVariables &primaryVars)
{
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
meanMolarMass_[phaseIdx] = 0;
}
// extract the mole fractions in both phases
Scalar sumFuga = 0;
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
fugacity_[compIdx] = primaryVars[fug0Idx + compIdx];
Valgrind::CheckDefined(fugacity_[compIdx]);
sumFuga += fugacity_[compIdx];
// convert the fugacities into mole fractions
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
moleFrac_[phaseIdx][compIdx] =
fugacity_[compIdx] /
FluidSystem::activityCoeff(phaseIdx,
compIdx,
temperature_,
phasePressure(phaseIdx),
*this);
Valgrind::CheckDefined(moleFrac_[phaseIdx][compIdx]);
// update the mean molar mass of each phase
meanMolarMass_[phaseIdx]
+= FluidSystem::molarMass(compIdx)*moleFrac_[phaseIdx][compIdx];
}
}
// make sure that the primary variables do not represent an
// unphysical state!
if (sumFuga < 1e-30) {
DUNE_THROW(NumericalProblem,
"Sum of component fugacities is too small: " << sumFuga);
}
// calculate the total concentration of each phase
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
// make sure that the mean molar mass within a phase is
// not "way off"
if (std::abs(meanMolarMass_[phaseIdx]) < 1e-30)
DUNE_THROW(NumericalProblem,
"Mean molar mass of phase "
<< phaseIdx
<< " is too small ("
<< meanMolarMass_[phaseIdx]
<< ")\n" << primaryVars);
// calculate the total concentration of the phase
phaseConcentration_[phaseIdx] =
FluidSystem::phaseDensity(phaseIdx,
temperature(),
phasePressure(phaseIdx),
*this)
/
meanMolarMass_[phaseIdx];
}
Valgrind::CheckDefined(fugacity_);
Valgrind::CheckDefined(moleFrac_);
Valgrind::CheckDefined(phaseConcentration_);
Valgrind::CheckDefined(meanMolarMass_);
Valgrind::CheckDefined(temperature_);
Valgrind::CheckDefined(phasePressure_);
Valgrind::CheckDefined(saturation_);
}
public:
/*!
* \brief Returns the saturation of a phase.
*/
Scalar saturation(int phaseIdx) const
{ return saturation_[phaseIdx]; }
/*!
* \brief Returns a vector of all phase saturations.
*/
const Scalar *saturations() const
{ return saturation_; }
/*!
* \brief Returns the molar fraction of a component in a fluid phase.
*/
Scalar moleFrac(int phaseIdx, int compIdx) const
{ return moleFrac_[phaseIdx][compIdx]; }
/*!
* \brief Returns the total concentration of a phase [mol / m^3].
*
* This is equivalent to the sum of all component concentrations.
*/
Scalar phaseConcentration(int phaseIdx) const
{ return phaseConcentration_[phaseIdx]; };
/*!
* \brief Returns the concentration of a component in a phase [mol / m^3].
*/
Scalar concentration(int phaseIdx, int compIdx) const
{ return moleFrac_[phaseIdx][compIdx]*phaseConcentration_[phaseIdx]; };
/*!
* \brief Returns the mass fraction of a component in a phase.
*/
Scalar massFrac(int phaseIdx, int compIdx) const
{
return
moleFrac_[phaseIdx][compIdx]*
FluidSystem::molarMass(compIdx)
/
meanMolarMass_[phaseIdx];
}
/*!
* \brief Returns the density of a phase [kg / m^3].
*/
Scalar density(int phaseIdx) const
{ return phaseConcentration_[phaseIdx]*meanMolarMass_[phaseIdx]; }
/*!
* \brief Returns mean molar mass of a phase [kg / mol].
*
* This is equivalent to the sum of all component molar masses
* weighted by their respective mole fraction.
*/
Scalar meanMolarMass(int phaseIdx) const
{ return meanMolarMass_[phaseIdx]; };
/*!
* \brief Returns the fugacity of a component [Pa].
*/
Scalar fugacity(int compIdx) const
{ return fugacity_[compIdx]; }
/*!
* \brief Returns the pressure of a fluid phase [Pa].
*/
Scalar phasePressure(int phaseIdx) const
{ return phasePressure_[phaseIdx]; }
/*!
* \brief Returns the capillary pressure [Pa]
*/
Scalar capillaryPressure(int phaseIdx) const
{ return phasePressure_[0] - phasePressure_[phaseIdx]; }
/*!
* \brief Returns the temperature of the fluids [K].
*
* Note that we assume thermodynamic equilibrium, so all fluids
* and the rock matrix exhibit the same temperature.
*/
Scalar temperature() const
{ return temperature_; };
public:
Scalar fugacity_[numComponents];
Scalar moleFrac_[numPhases][numComponents];
Scalar phaseConcentration_[numPhases];
Scalar meanMolarMass_[numPhases];
Scalar phasePressure_[numPhases];
Scalar saturation_[numPhases];
Scalar temperature_;
};
} // end namespace Opm
#endif // OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED

View File

@ -0,0 +1,11 @@
noinst_LTLIBRARIES = libblackoilfluid_noinst.la
libblackoilfluid_noinst_la_SOURCES = \
BlackoilPVT.cpp \
MiscibilityDead.cpp \
MiscibilityLiveGas.cpp \
MiscibilityLiveOil.cpp \
MiscibilityProps.cpp
include $(top_srcdir)/am/global-rules

View File

@ -0,0 +1,96 @@
//===========================================================================
//
// File: MiscibilityDead.cpp
//
// Created: Wed Feb 10 09:06:04 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
#include <algorithm>
#include "MiscibilityDead.hpp"
#include <dune/common/ErrorMacros.hpp>
#include <dune/common/linInt.hpp>
using namespace std;
using namespace Dune;
namespace samcode
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityDead::MiscibilityDead(const table_t& pvd_table)
: pvdx_(pvd_table)
{
const int region_number = 0;
if (pvd_table.size() != 1) {
THROW("More than one PVT-region");
}
// Convert units
const double bar = 1e5;
const double VISCOSITY_UNIT = 1e-3;
const int sz = pvdx_[region_number][0].size();
for (int i=0; i<sz; ++i) {
pvdx_[region_number][0][i] *= bar; // Pressure
pvdx_[region_number][2][i] *= VISCOSITY_UNIT;
}
// Interpolate 1/B
for (int i=0; i<sz; ++i) {
pvdx_[region_number][1][i] = 1.0/pvdx_[region_number][1][i];
}
}
// Destructor
MiscibilityDead::~MiscibilityDead()
{
}
double MiscibilityDead::getViscosity(int region, double press, const surfvol_t& /*surfvol*/) const
{
return linearInterpolationExtrap(pvdx_[region][0],
pvdx_[region][2], press);
}
double MiscibilityDead::B(int region, double press, const surfvol_t& /*surfvol*/) const
{
// Interpolate 1/B
return 1.0/linearInterpolationExtrap(pvdx_[region][0],
pvdx_[region][1], press);
//return linearInterpolationExtrap(pvdx_[region][0],
//pvdx_[region_number_][1], press);
}
double MiscibilityDead::dBdp(int region, double press, const surfvol_t& /*surfvol*/) const
{
// Interpolate 1/B
surfvol_t dummy_surfvol;
double Bg = B(region, press, dummy_surfvol);
return -Bg*Bg*
linearInterpolDerivative(pvdx_[region][0],
pvdx_[region][1], press);
//return linearInterpolDerivative(pvdx_[region][0],
// pvdx_[region][1], press);
}
double MiscibilityDead::R(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
{
return 0.0;
}
double MiscibilityDead::dRdp(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
{
return 0.0;
}
}

View File

@ -0,0 +1,47 @@
//===========================================================================
//
// File: MiscibilityDead.hpp
//
// Created: Wed Feb 10 09:05:47 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
#ifndef SINTEF_MISCIBILITYDEAD_HEADER
#define SINTEF_MISCIBILITYDEAD_HEADER
/** Class for immiscible dead oil and dry gas.
* Detailed description.
*/
#include "MiscibilityProps.hpp"
namespace samcode
{
class MiscibilityDead : public MiscibilityProps
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
MiscibilityDead(const table_t& pvd_table);
virtual ~MiscibilityDead();
virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const;
virtual double B(int region, double press, const surfvol_t& surfvol) const;
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const;
virtual double R(int region, double press, const surfvol_t& surfvol) const;
virtual double dRdp(int region, double press, const surfvol_t& surfvol) const;
private:
// PVT properties of dry gas or dead oil
table_t pvdx_;
};
}
#endif // SINTEF_MISCIBILITYDEAD_HEADER

View File

@ -0,0 +1,200 @@
//===========================================================================
//
// File: MiscibilityLiveGas.cpp
//
// Created: Wed Feb 10 09:21:53 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
#include <algorithm>
#include "MiscibilityLiveGas.hpp"
#include <dune/common/ErrorMacros.hpp>
#include <dune/common/linInt.hpp>
using namespace std;
using namespace Dune;
namespace samcode
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityLiveGas::MiscibilityLiveGas(const table_t& pvtg)
{
// GAS, PVTG
const double bar = 1e5;
const double VISCOSITY_UNIT = 1e-3;
const int region_number = 0;
if (pvtg.size() != 1) {
THROW("More than one PVD-region");
}
saturated_gas_table_.resize(4);
const int sz = pvtg[region_number].size();
for (int k=0; k<4; ++k) {
saturated_gas_table_[k].resize(sz);
}
for (int i=0; i<sz; ++i) {
saturated_gas_table_[0][i] = pvtg[region_number][i][0]*bar; // p
saturated_gas_table_[1][i] = pvtg[region_number][i][2]; // Bg
saturated_gas_table_[2][i] = pvtg[region_number][i][3]*VISCOSITY_UNIT; // mu_g
saturated_gas_table_[3][i] = pvtg[region_number][i][1]; // Rv
}
undersat_gas_tables_.resize(sz);
for (int i=0; i<sz; ++i) {
undersat_gas_tables_[i].resize(3);
int tsize = (pvtg[region_number][i].size() - 1)/3;
undersat_gas_tables_[i][0].resize(tsize);
undersat_gas_tables_[i][1].resize(tsize);
undersat_gas_tables_[i][2].resize(tsize);
for (int j=0, k=0; j<tsize; ++j) {
undersat_gas_tables_[i][0][j] = pvtg[region_number][i][++k]; // Rv
undersat_gas_tables_[i][1][j] = pvtg[region_number][i][++k]; // Bg
undersat_gas_tables_[i][2][j] = pvtg[region_number][i][++k]*VISCOSITY_UNIT; // mu_g
}
}
}
// Destructor
MiscibilityLiveGas::~MiscibilityLiveGas()
{
}
double MiscibilityLiveGas::getViscosity(int /*region*/, double press, const surfvol_t& surfvol) const
{
return miscible_gas(press, surfvol, 2, false);
}
// Vaporised oil-gas ratio.
double MiscibilityLiveGas::R(int /*region*/, double press, const surfvol_t& surfvol) const
{
if (surfvol[Liquid] == 0.0) {
return 0.0;
}
double R = linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[3], press);
double maxR = surfvol[Liquid]/surfvol[Vapour];
if (R < maxR ) { // Saturated case
return R;
} else {
return maxR; // Undersaturated case
}
}
// Vaporised oil-gas ratio derivative
double MiscibilityLiveGas::dRdp(int /*region*/, double press, const surfvol_t& surfvol) const
{
double R = linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[3], press);
double maxR = surfvol[Liquid]/surfvol[Vapour];
if (R < maxR ) { // Saturated case
return linearInterpolDerivative(saturated_gas_table_[0],
saturated_gas_table_[3],
press);
} else {
return 0.0; // Undersaturated case
}
}
double MiscibilityLiveGas::B(int /*region*/, double press, const surfvol_t& surfvol) const
{
if (surfvol[Vapour] == 0.0) return 1.0; // To handle no-gas case.
return miscible_gas(press, surfvol, 1, false);
}
double MiscibilityLiveGas::dBdp(int /*region*/, double press, const surfvol_t& surfvol) const
{
if (surfvol[Vapour] == 0.0) return 0.0; // To handle no-gas case.
return miscible_gas(press, surfvol, 1, true);
}
double MiscibilityLiveGas::miscible_gas(double press, const surfvol_t& surfvol, int item,
bool deriv) const
{
int section;
double R = linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[3], press,
section);
double maxR = surfvol[Liquid]/surfvol[Vapour];
if (deriv) {
if (R < maxR ) { // Saturated case
return linearInterpolDerivative(saturated_gas_table_[0],
saturated_gas_table_[item],
press);
} else { // Undersaturated case
int is = section;
if (undersat_gas_tables_[is][0].size() < 2) {
double val = (saturated_gas_table_[item][is+1]
- saturated_gas_table_[item][is]) /
(saturated_gas_table_[0][is+1] -
saturated_gas_table_[0][is]);
return val;
}
double val1 =
linearInterpolationExtrap(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
maxR);
double val2 =
linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
maxR);
double val = (val2 - val1)/
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
return val;
}
} else {
if (R < maxR ) { // Saturated case
return linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[item],
press);
} else { // Undersaturated case
int is = section;
// Extrapolate from first table section
if (is == 0 && press < saturated_gas_table_[0][0]) {
return linearInterpolationExtrap(undersat_gas_tables_[0][0],
undersat_gas_tables_[0][item],
maxR);
}
// Extrapolate from last table section
int ltp = saturated_gas_table_[0].size() - 1;
if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
return linearInterpolationExtrap(undersat_gas_tables_[ltp][0],
undersat_gas_tables_[ltp][item],
maxR);
}
// Interpolate between table sections
double w = (press - saturated_gas_table_[0][is]) /
(saturated_gas_table_[0][is+1] -
saturated_gas_table_[0][is]);
if (undersat_gas_tables_[is][0].size() < 2) {
double val = saturated_gas_table_[item][is] +
w*(saturated_gas_table_[item][is+1] -
saturated_gas_table_[item][is]);
return val;
}
double val1 =
linearInterpolationExtrap(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
maxR);
double val2 =
linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
maxR);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace samcode

View File

@ -0,0 +1,51 @@
//===========================================================================
//
// File: MiscibilityLiveGas.hpp
//
// Created: Wed Feb 10 09:21:26 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
#ifndef SINTEF_MISCIBILITYLIVEGAS_HEADER
#define SINTEF_MISCIBILITYLIVEGAS_HEADER
/** Class for miscible wet gas.
* Detailed description.
*/
#include "MiscibilityProps.hpp"
namespace samcode
{
class MiscibilityLiveGas : public MiscibilityProps
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
MiscibilityLiveGas(const table_t& pvto);
virtual ~MiscibilityLiveGas();
virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const;
virtual double R(int region, double press, const surfvol_t& surfvol) const;
virtual double dRdp(int region, double press, const surfvol_t& surfvol) const;
virtual double B(int region, double press, const surfvol_t& surfvol) const;
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const;
protected:
// item: 1=B 2=mu;
double miscible_gas(double press, const surfvol_t& surfvol, int item,
bool deriv = false) const;
// PVT properties of wet gas (with vaporised oil)
std::vector<std::vector<double> > saturated_gas_table_;
std::vector<std::vector<std::vector<double> > > undersat_gas_tables_;
};
}
#endif // SINTEF_MISCIBILITYLIVEGAS_HEADER

View File

@ -0,0 +1,203 @@
//===========================================================================
//
// File: MiscibiltyLiveOil.cpp
//
// Created: Wed Feb 10 09:08:25 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
#include <algorithm>
#include "MiscibilityLiveOil.hpp"
#include <dune/common/ErrorMacros.hpp>
#include <dune/common/linInt.hpp>
using namespace std;
using namespace Dune;
namespace samcode
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityLiveOil::MiscibilityLiveOil(const table_t& pvto)
{
// OIL, PVTO
const double bar = 1e5;
const double VISCOSITY_UNIT = 1e-3;
const int region_number = 0;
if (pvto.size() != 1) {
THROW("More than one PVD-region");
}
saturated_oil_table_.resize(4);
const int sz = pvto[region_number].size();
for (int k=0; k<4; ++k) {
saturated_oil_table_[k].resize(sz);
}
for (int i=0; i<sz; ++i) {
saturated_oil_table_[0][i] = pvto[region_number][i][1]*bar; // p
saturated_oil_table_[1][i] = 1.0/pvto[region_number][i][2]; // 1/Bo
saturated_oil_table_[2][i] = pvto[region_number][i][3]*VISCOSITY_UNIT; // mu_o
saturated_oil_table_[3][i] = pvto[region_number][i][0]; // Rs
}
undersat_oil_tables_.resize(sz);
for (int i=0; i<sz; ++i) {
undersat_oil_tables_[i].resize(3);
int tsize = (pvto[region_number][i].size() - 1)/3;
undersat_oil_tables_[i][0].resize(tsize);
undersat_oil_tables_[i][1].resize(tsize);
undersat_oil_tables_[i][2].resize(tsize);
for (int j=0, k=0; j<tsize; ++j) {
undersat_oil_tables_[i][0][j] = pvto[region_number][i][++k]*bar; // p
undersat_oil_tables_[i][1][j] = 1.0/pvto[region_number][i][++k]; // 1/Bo
undersat_oil_tables_[i][2][j] = pvto[region_number][i][++k]*VISCOSITY_UNIT; // mu_o
}
}
}
// Destructor
MiscibilityLiveOil::~MiscibilityLiveOil()
{
}
double MiscibilityLiveOil::getViscosity(int /*region*/, double press, const surfvol_t& surfvol) const
{
return miscible_oil(press, surfvol, 2, false);
}
// Dissolved gas-oil ratio
double MiscibilityLiveOil::R(int /*region*/, double press, const surfvol_t& surfvol) const
{
if (surfvol[Vapour] == 0.0) {
return 0.0;
}
double R = linearInterpolationExtrap(saturated_oil_table_[0],
saturated_oil_table_[3], press);
double maxR = surfvol[Vapour]/surfvol[Liquid];
if (R < maxR ) { // Saturated case
return R;
} else {
return maxR; // Undersaturated case
}
}
// Dissolved gas-oil ratio derivative
double MiscibilityLiveOil::dRdp(int /*region*/, double press, const surfvol_t& surfvol) const
{
double R = linearInterpolationExtrap(saturated_oil_table_[0],
saturated_oil_table_[3], press);
double maxR = surfvol[Vapour]/surfvol[Liquid];
if (R < maxR ) { // Saturated case
return linearInterpolDerivative(saturated_oil_table_[0],
saturated_oil_table_[3],
press);
} else {
return 0.0; // Undersaturated case
}
}
double MiscibilityLiveOil::B(int /*region*/, double press, const surfvol_t& surfvol) const
{
if (surfvol[Liquid] == 0.0) return 1.0; // To handle no-oil case.
return 1.0/miscible_oil(press, surfvol, 1, false);
}
double MiscibilityLiveOil::dBdp(int region, double press, const surfvol_t& surfvol) const
{
if (surfvol[Liquid] == 0.0) return 0.0; // To handle no-oil case.
double Bo = B(region, press, surfvol);
return -Bo*Bo*miscible_oil(press, surfvol, 1, true);
}
double MiscibilityLiveOil::miscible_oil(double press, const surfvol_t& surfvol,
int item, bool deriv) const
{
int section;
double R = linearInterpolationExtrap(saturated_oil_table_[0],
saturated_oil_table_[3],
press, section);
double maxR = surfvol[Vapour]/surfvol[Liquid];
if (deriv) {
if (R < maxR ) { // Saturated case
return linearInterpolDerivative(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[3], maxR);
if (undersat_oil_tables_[is][0].size() < 2) {
double val = (saturated_oil_table_[item][is+1]
- saturated_oil_table_[item][is]) /
(saturated_oil_table_[0][is+1] -
saturated_oil_table_[0][is]);
return val;
}
double w = (maxR - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
double val1 =
linearInterpolDerivative(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolDerivative(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
} else {
if (R < maxR ) { // Saturated case
return linearInterpolationExtrap(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[3], maxR);
// Extrapolate from first table section
if (is == 0 && press < saturated_oil_table_[0][0]) {
return linearInterpolationExtrap(undersat_oil_tables_[0][0],
undersat_oil_tables_[0][item],
maxR);
}
// Extrapolate from last table section
int ltp = saturated_oil_table_[0].size() - 1;
if (is+1 == ltp && press > saturated_oil_table_[0][ltp]) {
return linearInterpolationExtrap(undersat_oil_tables_[ltp][0],
undersat_oil_tables_[ltp][item],
maxR);
}
// Interpolate between table sections
double w = (maxR - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] -
saturated_oil_table_[3][is]);
if (undersat_oil_tables_[is][0].size() < 2) {
double val = saturated_oil_table_[item][is] +
w*(saturated_oil_table_[item][is+1] -
saturated_oil_table_[item][is]);
return val;
}
double val1 =
linearInterpolationExtrap(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolationExtrap(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace samcode

View File

@ -0,0 +1,51 @@
//===========================================================================
//
// File: MiscibilityLiveOil.hpp
//
// Created: Wed Feb 10 09:08:09 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
#ifndef SINTEF_MISCIBILITYLIVEOIL_HEADER
#define SINTEF_MISCIBILITYLIVEOIL_HEADER
/** Class for miscible live oil.
* Detailed description.
*/
#include "MiscibilityProps.hpp"
namespace samcode
{
class MiscibilityLiveOil : public MiscibilityProps
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
MiscibilityLiveOil(const table_t& pvto);
virtual ~MiscibilityLiveOil();
virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const;
virtual double R (int region, double press, const surfvol_t& surfvol) const;
virtual double dRdp(int region, double press, const surfvol_t& surfvol) const;
virtual double B (int region, double press, const surfvol_t& surfvol) const;
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const;
protected:
// item: 1=B 2=mu;
double miscible_oil(double press, const surfvol_t& surfvol, int item,
bool deriv = false) const;
// PVT properties of live oil (with dissolved gas)
std::vector<std::vector<double> > saturated_oil_table_;
std::vector<std::vector<std::vector<double> > > undersat_oil_tables_;
};
}
#endif // SINTEF_MISCIBILITYLIVEOIL_HEADER

View File

@ -0,0 +1,34 @@
//===========================================================================
//
// File: MiscibilityProps.cpp
//
// Created: Wed Feb 10 09:05:05 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
#include "MiscibilityProps.hpp"
using namespace std;
namespace samcode
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityProps::MiscibilityProps()
{
}
MiscibilityProps::~MiscibilityProps()
{
}
} // namespace samcode

View File

@ -0,0 +1,46 @@
//===========================================================================
//
// File: MiscibilityProps.hpp
//
// Created: Wed Feb 10 09:04:35 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
#ifndef SINTEF_MISCIBILITYPROPS_HEADER
#define SINTEF_MISCIBILITYPROPS_HEADER
/** Base class for properties of fluids and rocks.
* Detailed description.
*/
#include <vector>
#include <tr1/array>
namespace samcode
{
class MiscibilityProps
{
public:
typedef std::tr1::array<double, 3> surfvol_t;
enum PhaseNames { Aqua = 0, Liquid = 1, Vapour = 2 };
MiscibilityProps();
virtual ~MiscibilityProps();
virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const = 0;
virtual double B (int region, double press, const surfvol_t& surfvol) const = 0;
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const = 0;
virtual double R (int region, double press, const surfvol_t& surfvol) const = 0;
virtual double dRdp(int region, double press, const surfvol_t& surfvol) const = 0;
};
} // namespace samcode
#endif // SINTEF_MISCIBILITYPROPS_HEADER

View File

@ -0,0 +1,81 @@
//===========================================================================
//
// File: MiscibilityWater.hpp
//
// Created: Tue May 18 10:26:13 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Bjørn Spjelkavik <bsp@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
#ifndef OPENRS_MISCIBILITYWATER_HEADER
#define OPENRS_MISCIBILITYWATER_HEADER
#include "MiscibilityProps.hpp"
#include <dune/common/ErrorMacros.hpp>
// Forward declaration.
class PVTW;
namespace samcode
{
class MiscibilityWater : public MiscibilityProps
{
public:
typedef std::vector<std::vector<double> > table_t;
MiscibilityWater(const table_t& pvtw)
{
const int region_number = 0;
if (pvtw.size() != 1) {
THROW("More than one PVD-region");
}
double b = pvtw[region_number][1];
double comp = pvtw[region_number][2];
double visc = pvtw[region_number][3];
const double VISCOSITY_UNIT = 1e-3;
if (b == 1.0 && comp == 0) {
viscosity_ = visc*VISCOSITY_UNIT;
} else {
THROW("Compressible water not implemented.");
}
}
MiscibilityWater(double visc)
: viscosity_(visc)
{
}
virtual ~MiscibilityWater()
{
}
virtual double getViscosity(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
{
return viscosity_;
}
virtual double B(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
{
return 1.0;
}
virtual double dBdp(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
{
return 0.0;
}
virtual double R(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
{
return 0.0;
}
virtual double dRdp(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
{
return 0.0;
}
private:
double viscosity_;
};
}
#endif // OPENRS_MISCIBILITYWATER_HEADER