36#ifndef OPENRS_RESERVOIRPROPERTYCAPILLARY_IMPL_HEADER
37#define OPENRS_RESERVOIRPROPERTYCAPILLARY_IMPL_HEADER
47 return relPermFirstPhase(cell_index, saturation) / Super::viscosity1_;
54 return relPermSecondPhase(cell_index, saturation) / Super::viscosity2_;
62 double& phase_mob)
const
64 if (phase_index == 0) {
65 phase_mob = mobilityFirstPhase(cell_index, saturation);
67 assert(phase_index == 1);
68 phase_mob = mobilitySecondPhase(cell_index, saturation);
76 double l1 = mobilityFirstPhase(cell_index, saturation);
77 double l2 = mobilitySecondPhase(cell_index, saturation);
85 double l1 = mobilityFirstPhase(cell_index, saturation);
86 double l2 = mobilitySecondPhase(cell_index, saturation);
92 template<
class Vector>
96 mobility[0] = mobilityFirstPhase(cell_index, saturation);
97 mobility[1] = mobilitySecondPhase(cell_index, saturation);
102 template <
class Vector>
105 Vector& dmob)
const {
107 dmob[0] = relPermFirstPhaseDeriv (c, s) / Super::viscosity1_;
108 dmob[3] = - relPermSecondPhaseDeriv(c, s) / Super::viscosity2_;
109 dmob[1] = dmob[2] = 0;
120 if (Super::rock_.size() > 0) {
121 const int region = Super::cell_to_rock_[cell_index];
122 assert (region <
int(Super::rock_.size()));
124 Super::rock_[region].krw(saturation, res);
129 return saturation * saturation;
135 ReservoirPropertyCapillary<dim>::
136 relPermFirstPhaseDeriv(
int cell_index,
double saturation)
const
138 if (Super::rock_.size() > 0) {
139 const int region = Super::cell_to_rock_[cell_index];
140 assert (region <
int(Super::rock_.size()));
142 Super::rock_[region].dkrw(saturation, res);
147 return 2 * saturation;
155 double ReservoirPropertyCapillary<dim>::relPermSecondPhase(
int cell_index,
double saturation)
const
157 if (Super::rock_.size() > 0) {
158 const int region = Super::cell_to_rock_[cell_index];
159 assert (region <
int(Super::rock_.size()));
161 Super::rock_[region].kro(saturation, res);
166 return (1 - saturation) * (1 - saturation);
172 ReservoirPropertyCapillary<dim>::
173 relPermSecondPhaseDeriv(
int cell_index,
double saturation)
const
175 if (Super::rock_.size() > 0) {
176 const int region = Super::cell_to_rock_[cell_index];
177 assert (region <
int(Super::rock_.size()));
179 Super::rock_[region].dkro(saturation, res);
184 return - 2 * (1 - saturation);
191 void ReservoirPropertyCapillary<dim>::cflFracFlows(
int rock,
double s,
double& ff_first,
double& ff_gravity)
const
195 const int cell_index = 0;
196 double l1 = mobilityFirstPhase(cell_index, s);
197 double l2 = mobilitySecondPhase(cell_index, s);
198 ff_first = l1/(l1 + l2);
199 ff_gravity = l1*l2/(l1 + l2);
202 Super::rock_[rock].krw(s, krw);
203 Super::rock_[rock].kro(s, kro);
204 double l1 = krw/Super::viscosity1_;
205 double l2 = kro/Super::viscosity2_;
206 ff_first = l1/(l1 + l2);
207 ff_gravity = l1*l2/(l1 + l2);
215 std::array<double, 3> ReservoirPropertyCapillary<dim>::computeSingleRockCflFactors(
int rock,
double min_perm,
double max_poro)
const
218 OwnCMatrix min_perm_matrix(dim, dim, (
double*)0);
219 eye(min_perm_matrix);
220 min_perm_matrix *= min_perm;
224 double delta = 1.0/double(N - 1);
225 double last_ff1, last_ffg;
226 double max_der1 = -1e100;
227 double max_derg = -1e100;
228 cflFracFlows(rock, 0.0, last_ff1, last_ffg);
229 double max_ffg = last_ffg;
230 double max_derpc = rock == -1 ? 0.0 :
231 std::fabs(Super::rock_[rock].capPressDeriv(min_perm_matrix, max_poro, 0.0));
232 for (
int i = 1; i < N; ++i) {
233 double s = double(i)*delta;
235 cflFracFlows(rock, s, ff1, ffg);
236 double est_deriv_ff1 = std::fabs(ff1 - last_ff1)/delta;
237 double est_deriv_ffg = std::fabs(ffg - last_ffg)/delta;
238 max_der1 = std::max(max_der1, est_deriv_ff1);
239 max_derg = std::max(max_derg, est_deriv_ffg);
240 max_ffg = std::max(max_ffg, ffg);
241 max_derpc = rock == -1 ? 0.0 :
242 std::max(max_derpc, std::fabs(Super::rock_[rock].capPressDeriv(min_perm_matrix, max_poro, s)));
246 std::array<double, 3> retval = {{ 1.0/max_der1, 1.0/max_derg, max_ffg*max_derpc }};
256 if (Super::rock_.empty()) {
257 std::array<double, 3> fac = computeSingleRockCflFactors(-1, 0.0, 0.0);
258 Super::cfl_factor_ = fac[0];
259 Super::cfl_factor_gravity_ = fac[1];
260 Super::cfl_factor_capillary_ = fac[2];
263 std::vector<double> min_perm(Super::rock_.size(), 1e100);
264 std::vector<double> max_poro(Super::rock_.size(), 0.0);
265 int num_cells = Super::porosity_.size();
266 for (
int c = 0; c < num_cells; ++c) {
267 int r = Super::cell_to_rock_[c];
268 min_perm[r] = std::min(min_perm[r],
trace(Super::permeability(c))/
double(dim));
269 max_poro[r] = std::max(max_poro[r], Super::porosity(c));
271 Super::cfl_factor_ = 1e100;
272 Super::cfl_factor_gravity_ = 1e100;
273 Super::cfl_factor_capillary_ = 0.0;
274 for (
int r = 0; r < int(Super::rock_.size()); ++r) {
275 std::array<double, 3> fac = computeSingleRockCflFactors(r, min_perm[r], max_poro[r]);
276 Super::cfl_factor_ = std::min(Super::cfl_factor_, fac[0]);
277 Super::cfl_factor_gravity_ = std::min(Super::cfl_factor_gravity_, fac[1]);
278 Super::cfl_factor_capillary_ = std::max(Super::cfl_factor_capillary_, fac[2]);
A property class for incompressible two-phase flow.
Definition: ReservoirPropertyCapillary.hpp:80
void computeCflFactors()
Computes cfl factors. Called from ReservoirPropertyCommon::init().
Definition: ReservoirPropertyCapillary_impl.hpp:254
void phaseMobilities(int cell_index, double saturation, Vector &mobility) const
Mobilities for both phases.
Definition: ReservoirPropertyCapillary_impl.hpp:93
double mobilitySecondPhase(int cell_index, double saturation) const
Mobility of second (oil) phase.
Definition: ReservoirPropertyCapillary_impl.hpp:52
void phaseMobilitiesDeriv(int c, double s, Vector &dmob) const
Definition: ReservoirPropertyCapillary_impl.hpp:104
void phaseMobility(int phase_index, int cell_index, double saturation, double &phase_mob) const
Phase mobility.
Definition: ReservoirPropertyCapillary_impl.hpp:59
double mobilityFirstPhase(int cell_index, double saturation) const
Mobility of first (water) phase.
Definition: ReservoirPropertyCapillary_impl.hpp:45
double totalMobility(int cell_index, double saturation) const
Total mobility.
Definition: ReservoirPropertyCapillary_impl.hpp:74
double fractionalFlow(int cell_index, double saturation) const
Fractional flow (of the first phase).
Definition: ReservoirPropertyCapillary_impl.hpp:83
Definition: BlackoilFluid.hpp:32
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition: Matrix.hpp:639
void eye(Matrix &A)
Make an identity.
Definition: Matrix.hpp:618
FullMatrix< double, OwnData, COrdering > OwnCMatrix
Convenience typedefs for C-ordered.
Definition: Matrix.hpp:580