Loading...
Searching...
No Matches
ReservoirPropertyCapillary_impl.hpp
//===========================================================================
//
// File: ReservoirPropertyCapillary_impl.hpp
//
// Created: Thu Oct 22 20:16:15 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS 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.
OpenRS 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 OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_RESERVOIRPROPERTYCAPILLARY_IMPL_HEADER
#define OPENRS_RESERVOIRPROPERTYCAPILLARY_IMPL_HEADER
namespace Opm
{
template <int dim>
double ReservoirPropertyCapillary<dim>::mobilityFirstPhase(int cell_index, double saturation) const
{
return relPermFirstPhase(cell_index, saturation) / Super::viscosity1_;
}
template <int dim>
double ReservoirPropertyCapillary<dim>::mobilitySecondPhase(int cell_index, double saturation) const
{
return relPermSecondPhase(cell_index, saturation) / Super::viscosity2_;
}
template <int dim>
void ReservoirPropertyCapillary<dim>::phaseMobility(int phase_index,
int cell_index,
double saturation,
double& phase_mob) const
{
if (phase_index == 0) {
phase_mob = mobilityFirstPhase(cell_index, saturation);
} else {
assert(phase_index == 1);
phase_mob = mobilitySecondPhase(cell_index, saturation);
}
}
template <int dim>
double ReservoirPropertyCapillary<dim>::totalMobility(int cell_index, double saturation) const
{
double l1 = mobilityFirstPhase(cell_index, saturation);
double l2 = mobilitySecondPhase(cell_index, saturation);
return l1 + l2;
}
template <int dim>
double ReservoirPropertyCapillary<dim>::fractionalFlow(int cell_index, double saturation) const
{
double l1 = mobilityFirstPhase(cell_index, saturation);
double l2 = mobilitySecondPhase(cell_index, saturation);
return l1/(l1 + l2);
}
template <int dim>
template<class Vector>
void ReservoirPropertyCapillary<dim>::phaseMobilities(int cell_index, double saturation, Vector& mobility) const
{
//assert (mobility.size() >= Super::NumberOfPhases);
mobility[0] = mobilityFirstPhase(cell_index, saturation);
mobility[1] = mobilitySecondPhase(cell_index, saturation);
}
template <int dim>
template <class Vector>
void
ReservoirPropertyCapillary<dim>::phaseMobilitiesDeriv(int c, double s,
Vector& dmob) const {
dmob[0] = relPermFirstPhaseDeriv (c, s) / Super::viscosity1_;
dmob[3] = - relPermSecondPhaseDeriv(c, s) / Super::viscosity2_;
dmob[1] = dmob[2] = 0;
}
// ------ Private methods ------
template <int dim>
double ReservoirPropertyCapillary<dim>::relPermFirstPhase(int cell_index, double saturation) const
{
if (Super::rock_.size() > 0) {
const int region = Super::cell_to_rock_[cell_index];
assert (region < int(Super::rock_.size()));
double res;
Super::rock_[region].krw(saturation, res);
return res;
} else {
// HACK ALERT!
// Use quadratic rel-perm if no known rock table exists.
return saturation * saturation;
}
}
template <int dim>
double
ReservoirPropertyCapillary<dim>::
relPermFirstPhaseDeriv(int cell_index, double saturation) const
{
if (Super::rock_.size() > 0) {
const int region = Super::cell_to_rock_[cell_index];
assert (region < int(Super::rock_.size()));
double res;
Super::rock_[region].dkrw(saturation, res);
return res;
} else {
// HACK ALERT!
// Use quadratic rel-perm if no known rock table exists.
return 2 * saturation;
}
}
template <int dim>
double ReservoirPropertyCapillary<dim>::relPermSecondPhase(int cell_index, double saturation) const
{
if (Super::rock_.size() > 0) {
const int region = Super::cell_to_rock_[cell_index];
assert (region < int(Super::rock_.size()));
double res;
Super::rock_[region].kro(saturation, res);
return res;
} else {
// HACK ALERT!
// Use quadratic rel-perm if no known rock table exists.
return (1 - saturation) * (1 - saturation);
}
}
template <int dim>
double
ReservoirPropertyCapillary<dim>::
relPermSecondPhaseDeriv(int cell_index, double saturation) const
{
if (Super::rock_.size() > 0) {
const int region = Super::cell_to_rock_[cell_index];
assert (region < int(Super::rock_.size()));
double res;
Super::rock_[region].dkro(saturation, res);
return res;
} else {
// HACK ALERT!
// Use quadratic rel-perm if no known rock table exists.
return - 2 * (1 - saturation);
}
}
template <int dim>
void ReservoirPropertyCapillary<dim>::cflFracFlows(int rock, double s, double& ff_first, double& ff_gravity) const
{
if (rock == -1) {
// No rock dependency, we might just as well use the first cell.
const int cell_index = 0;
double l1 = mobilityFirstPhase(cell_index, s);
double l2 = mobilitySecondPhase(cell_index, s);
ff_first = l1/(l1 + l2);
ff_gravity = l1*l2/(l1 + l2);
} else {
double krw, kro;
Super::rock_[rock].krw(s, krw);
Super::rock_[rock].kro(s, kro);
double l1 = krw/Super::viscosity1_;
double l2 = kro/Super::viscosity2_;
ff_first = l1/(l1 + l2);
ff_gravity = l1*l2/(l1 + l2);
}
}
template <int dim>
std::array<double, 3> ReservoirPropertyCapillary<dim>::computeSingleRockCflFactors(int rock, double min_perm, double max_poro) const
{
// Make min_perm matrix.
OwnCMatrix min_perm_matrix(dim, dim, (double*)0);
eye(min_perm_matrix);
min_perm_matrix *= min_perm;
// Sample values at many saturation points.
const int N = 257;
double delta = 1.0/double(N - 1);
double last_ff1, last_ffg;
double max_der1 = -1e100;
double max_derg = -1e100;
cflFracFlows(rock, 0.0, last_ff1, last_ffg);
double max_ffg = last_ffg;
double max_derpc = rock == -1 ? 0.0 :
std::fabs(Super::rock_[rock].capPressDeriv(min_perm_matrix, max_poro, 0.0));
for (int i = 1; i < N; ++i) {
double s = double(i)*delta;
double ff1, ffg;
cflFracFlows(rock, s, ff1, ffg);
double est_deriv_ff1 = std::fabs(ff1 - last_ff1)/delta;
double est_deriv_ffg = std::fabs(ffg - last_ffg)/delta;
max_der1 = std::max(max_der1, est_deriv_ff1);
max_derg = std::max(max_derg, est_deriv_ffg);
max_ffg = std::max(max_ffg, ffg);
max_derpc = rock == -1 ? 0.0 :
std::max(max_derpc, std::fabs(Super::rock_[rock].capPressDeriv(min_perm_matrix, max_poro, s)));
last_ff1 = ff1;
last_ffg = ffg;
}
std::array<double, 3> retval = {{ 1.0/max_der1, 1.0/max_derg, max_ffg*max_derpc }};
return retval;
}
template <int dim>
void ReservoirPropertyCapillary<dim>::computeCflFactors()
{
if (Super::rock_.empty()) {
std::array<double, 3> fac = computeSingleRockCflFactors(-1, 0.0, 0.0);
Super::cfl_factor_ = fac[0];
Super::cfl_factor_gravity_ = fac[1];
Super::cfl_factor_capillary_ = fac[2];
} else {
// Compute min perm and max poro per rock (for J-scaling cap pressure funcs).
std::vector<double> min_perm(Super::rock_.size(), 1e100);
std::vector<double> max_poro(Super::rock_.size(), 0.0);
int num_cells = Super::porosity_.size();
for (int c = 0; c < num_cells; ++c) {
int r = Super::cell_to_rock_[c];
min_perm[r] = std::min(min_perm[r], trace(Super::permeability(c))/double(dim));
max_poro[r] = std::max(max_poro[r], Super::porosity(c));
}
Super::cfl_factor_ = 1e100;
Super::cfl_factor_gravity_ = 1e100;
Super::cfl_factor_capillary_ = 0.0;
for (int r = 0; r < int(Super::rock_.size()); ++r) {
std::array<double, 3> fac = computeSingleRockCflFactors(r, min_perm[r], max_poro[r]);
Super::cfl_factor_ = std::min(Super::cfl_factor_, fac[0]);
Super::cfl_factor_gravity_ = std::min(Super::cfl_factor_gravity_, fac[1]);
Super::cfl_factor_capillary_ = std::max(Super::cfl_factor_capillary_, fac[2]);
}
}
}
} // namespace Opm
#endif // OPENRS_RESERVOIRPROPERTYCAPILLARY_IMPL_HEADER
Opm::ReservoirPropertyCapillary
A property class for incompressible two-phase flow.
Definition ReservoirPropertyCapillary.hpp:80
Opm::ReservoirPropertyCapillary::computeCflFactors
void computeCflFactors()
Computes cfl factors. Called from ReservoirPropertyCommon::init().
Definition ReservoirPropertyCapillary_impl.hpp:254
Opm::ReservoirPropertyCapillary::phaseMobilities
void phaseMobilities(int cell_index, double saturation, Vector &mobility) const
Mobilities for both phases.
Definition ReservoirPropertyCapillary_impl.hpp:93
Opm::ReservoirPropertyCapillary::mobilitySecondPhase
double mobilitySecondPhase(int cell_index, double saturation) const
Mobility of second (oil) phase.
Definition ReservoirPropertyCapillary_impl.hpp:52
Opm::ReservoirPropertyCapillary::phaseMobility
void phaseMobility(int phase_index, int cell_index, double saturation, double &phase_mob) const
Phase mobility.
Definition ReservoirPropertyCapillary_impl.hpp:59
Opm::ReservoirPropertyCapillary::mobilityFirstPhase
double mobilityFirstPhase(int cell_index, double saturation) const
Mobility of first (water) phase.
Definition ReservoirPropertyCapillary_impl.hpp:45
Opm::ReservoirPropertyCapillary::totalMobility
double totalMobility(int cell_index, double saturation) const
Total mobility.
Definition ReservoirPropertyCapillary_impl.hpp:74
Opm::ReservoirPropertyCapillary::fractionalFlow
double fractionalFlow(int cell_index, double saturation) const
Fractional flow (of the first phase).
Definition ReservoirPropertyCapillary_impl.hpp:83
Opm::ReservoirPropertyCommon< dim, ReservoirPropertyCapillary< dim >, RockJfunc >::porosity
double porosity(int cell_index) const
Opm::ReservoirPropertyCommon< dim, ReservoirPropertyCapillary< dim >, RockJfunc >::permeability
PermTensor permeability(int cell_index) const
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::trace
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition Matrix.hpp:638
Opm::eye
void eye(Matrix &A)
Make an identity.
Definition Matrix.hpp:617
Opm::OwnCMatrix
FullMatrix< double, OwnData, COrdering > OwnCMatrix
Convenience typedefs for C-ordered.
Definition Matrix.hpp:579