ReservoirPropertyCapillary_impl.hpp
Go to the documentation of this file.
1//===========================================================================
2//
3// File: ReservoirPropertyCapillary_impl.hpp
4//
5// Created: Thu Oct 22 20:16:15 2009
6//
7// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8// Bård Skaflestad <bard.skaflestad@sintef.no>
9//
10// $Date$
11//
12// $Revision$
13//
14//===========================================================================
15
16/*
17 Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
18 Copyright 2009, 2010 Statoil ASA.
19
20 This file is part of The Open Reservoir Simulator Project (OpenRS).
21
22 OpenRS is free software: you can redistribute it and/or modify
23 it under the terms of the GNU General Public License as published by
24 the Free Software Foundation, either version 3 of the License, or
25 (at your option) any later version.
26
27 OpenRS is distributed in the hope that it will be useful,
28 but WITHOUT ANY WARRANTY; without even the implied warranty of
29 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 GNU General Public License for more details.
31
32 You should have received a copy of the GNU General Public License
33 along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
34*/
35
36#ifndef OPENRS_RESERVOIRPROPERTYCAPILLARY_IMPL_HEADER
37#define OPENRS_RESERVOIRPROPERTYCAPILLARY_IMPL_HEADER
38
39
40namespace Opm
41{
42
43
44 template <int dim>
45 double ReservoirPropertyCapillary<dim>::mobilityFirstPhase(int cell_index, double saturation) const
46 {
47 return relPermFirstPhase(cell_index, saturation) / Super::viscosity1_;
48 }
49
50
51 template <int dim>
52 double ReservoirPropertyCapillary<dim>::mobilitySecondPhase(int cell_index, double saturation) const
53 {
54 return relPermSecondPhase(cell_index, saturation) / Super::viscosity2_;
55 }
56
57
58 template <int dim>
59 void ReservoirPropertyCapillary<dim>::phaseMobility(int phase_index,
60 int cell_index,
61 double saturation,
62 double& phase_mob) const
63 {
64 if (phase_index == 0) {
65 phase_mob = mobilityFirstPhase(cell_index, saturation);
66 } else {
67 assert(phase_index == 1);
68 phase_mob = mobilitySecondPhase(cell_index, saturation);
69 }
70 }
71
72
73 template <int dim>
74 double ReservoirPropertyCapillary<dim>::totalMobility(int cell_index, double saturation) const
75 {
76 double l1 = mobilityFirstPhase(cell_index, saturation);
77 double l2 = mobilitySecondPhase(cell_index, saturation);
78 return l1 + l2;
79 }
80
81
82 template <int dim>
83 double ReservoirPropertyCapillary<dim>::fractionalFlow(int cell_index, double saturation) const
84 {
85 double l1 = mobilityFirstPhase(cell_index, saturation);
86 double l2 = mobilitySecondPhase(cell_index, saturation);
87 return l1/(l1 + l2);
88 }
89
90
91 template <int dim>
92 template<class Vector>
93 void ReservoirPropertyCapillary<dim>::phaseMobilities(int cell_index, double saturation, Vector& mobility) const
94 {
95 //assert (mobility.size() >= Super::NumberOfPhases);
96 mobility[0] = mobilityFirstPhase(cell_index, saturation);
97 mobility[1] = mobilitySecondPhase(cell_index, saturation);
98 }
99
100
101 template <int dim>
102 template <class Vector>
103 void
104 ReservoirPropertyCapillary<dim>::phaseMobilitiesDeriv(int c, double s,
105 Vector& dmob) const {
106
107 dmob[0] = relPermFirstPhaseDeriv (c, s) / Super::viscosity1_;
108 dmob[3] = - relPermSecondPhaseDeriv(c, s) / Super::viscosity2_;
109 dmob[1] = dmob[2] = 0;
110 }
111
112
113
114 // ------ Private methods ------
115
116
117 template <int dim>
118 double ReservoirPropertyCapillary<dim>::relPermFirstPhase(int cell_index, double saturation) const
119 {
120 if (Super::rock_.size() > 0) {
121 const int region = Super::cell_to_rock_[cell_index];
122 assert (region < int(Super::rock_.size()));
123 double res;
124 Super::rock_[region].krw(saturation, res);
125 return res;
126 } else {
127 // HACK ALERT!
128 // Use quadratic rel-perm if no known rock table exists.
129 return saturation * saturation;
130 }
131 }
132
133 template <int dim>
134 double
135 ReservoirPropertyCapillary<dim>::
136 relPermFirstPhaseDeriv(int cell_index, double saturation) const
137 {
138 if (Super::rock_.size() > 0) {
139 const int region = Super::cell_to_rock_[cell_index];
140 assert (region < int(Super::rock_.size()));
141 double res;
142 Super::rock_[region].dkrw(saturation, res);
143 return res;
144 } else {
145 // HACK ALERT!
146 // Use quadratic rel-perm if no known rock table exists.
147 return 2 * saturation;
148 }
149 }
150
151
152
153
154 template <int dim>
155 double ReservoirPropertyCapillary<dim>::relPermSecondPhase(int cell_index, double saturation) const
156 {
157 if (Super::rock_.size() > 0) {
158 const int region = Super::cell_to_rock_[cell_index];
159 assert (region < int(Super::rock_.size()));
160 double res;
161 Super::rock_[region].kro(saturation, res);
162 return res;
163 } else {
164 // HACK ALERT!
165 // Use quadratic rel-perm if no known rock table exists.
166 return (1 - saturation) * (1 - saturation);
167 }
168 }
169
170 template <int dim>
171 double
172 ReservoirPropertyCapillary<dim>::
173 relPermSecondPhaseDeriv(int cell_index, double saturation) const
174 {
175 if (Super::rock_.size() > 0) {
176 const int region = Super::cell_to_rock_[cell_index];
177 assert (region < int(Super::rock_.size()));
178 double res;
179 Super::rock_[region].dkro(saturation, res);
180 return res;
181 } else {
182 // HACK ALERT!
183 // Use quadratic rel-perm if no known rock table exists.
184 return - 2 * (1 - saturation);
185 }
186 }
187
188
189
190 template <int dim>
191 void ReservoirPropertyCapillary<dim>::cflFracFlows(int rock, double s, double& ff_first, double& ff_gravity) const
192 {
193 if (rock == -1) {
194 // No rock dependency, we might just as well use the first cell.
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);
200 } else {
201 double krw, kro;
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);
208 }
209 }
210
211
212
213
214 template <int dim>
215 std::array<double, 3> ReservoirPropertyCapillary<dim>::computeSingleRockCflFactors(int rock, double min_perm, double max_poro) const
216 {
217 // Make min_perm matrix.
218 OwnCMatrix min_perm_matrix(dim, dim, (double*)0);
219 eye(min_perm_matrix);
220 min_perm_matrix *= min_perm;
221
222 // Sample values at many saturation points.
223 const int N = 257;
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;
234 double ff1, ffg;
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)));
243 last_ff1 = ff1;
244 last_ffg = ffg;
245 }
246 std::array<double, 3> retval = {{ 1.0/max_der1, 1.0/max_derg, max_ffg*max_derpc }};
247 return retval;
248 }
249
250
251
252
253 template <int dim>
254 void ReservoirPropertyCapillary<dim>::computeCflFactors()
255 {
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];
261 } else {
262 // Compute min perm and max poro per rock (for J-scaling cap pressure funcs).
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));
270 }
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]);
279 }
280 }
281 }
282
283
284
285
286
287} // namespace Opm
288
289
290#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::phaseMobilitiesDeriv
void phaseMobilitiesDeriv(int c, double s, Vector &dmob) const
Definition: ReservoirPropertyCapillary_impl.hpp:104
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
Definition: BlackoilFluid.hpp:32
Opm::trace
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition: Matrix.hpp:639
Opm::eye
void eye(Matrix &A)
Make an identity.
Definition: Matrix.hpp:618
Opm::OwnCMatrix
FullMatrix< double, OwnData, COrdering > OwnCMatrix
Convenience typedefs for C-ordered.
Definition: Matrix.hpp:580