20#ifndef OPM_BLACKOILCO2PVT_HEADER_INCLUDED
21#define OPM_BLACKOILCO2PVT_HEADER_INCLUDED
22#define OPM_BLACKOILPVT_HEADER_INCLUDED
24#define OPM_DEPRECATED __attribute__((deprecated))
25#define OPM_DEPRECATED_MSG(msg) __attribute__((deprecated))
27#include <opm/material/fluidsystems/BrineCO2FluidSystem.hpp>
28#include <opm/material/fluidstates/CompositionalFluidState.hpp>
29#include <opm/material/constraintsolvers/MiscibleMultiPhaseComposition.hpp>
30#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
32#include <opm/common/Exceptions.hpp>
33#include <opm/common/ErrorMacros.hpp>
38#include <opm/parser/eclipse/Deck/Deck.hpp>
51 typedef Opm::FluidSystems::BrineCO2<double, Opm::Benchmark3::CO2Tables>
FluidSystem;
54 void init(Opm::DeckConstPtr deck);
65 double B (
double press,
68 double dBdp(
double press,
71 double R (
double press,
74 double dRdp(
double press,
78 void getViscosity(
const std::vector<PhaseVec>& pressures,
79 const std::vector<CompVec>& surfvol,
80 std::vector<PhaseVec>& output)
const;
81 void B(
const std::vector<PhaseVec>& pressures,
82 const std::vector<CompVec>& surfvol,
83 std::vector<PhaseVec>& output)
const;
84 void dBdp(
const std::vector<PhaseVec>& pressures,
85 const std::vector<CompVec>& surfvol,
86 std::vector<PhaseVec>& output_B,
87 std::vector<PhaseVec>& output_dBdp)
const;
88 void R(
const std::vector<PhaseVec>& pressures,
89 const std::vector<CompVec>& surfvol,
90 std::vector<PhaseVec>& output)
const;
91 void dRdp(
const std::vector<PhaseVec>& pressures,
92 const std::vector<CompVec>& surfvol,
93 std::vector<PhaseVec>& output_R,
94 std::vector<PhaseVec>& output_dRdp)
const;
102 Dune::FieldVector<double,2> density;
103 Dune::FieldMatrix<double,2,2> massfrac;
104 Dune::FieldVector<double,2> phaseVolume;
105 Dune::FieldVector<double,2> phaseViscosity;
108 void computeState(SubState& ss,
double zBrine,
double zCO2,
double pressure)
const;
110 wPhase = FluidSystem::liquidPhaseIdx,
111 nPhase = FluidSystem::gasPhaseIdx,
113 wComp = FluidSystem::BrineIdx,
114 nComp = FluidSystem::CO2Idx
123 surfaceDensities_[
Water] = 1000.;
124 surfaceDensities_[
Gas] = 2.0;
125 surfaceDensities_[
Oil] = 1000.;
134 std::cout <<
"\n Generating pvt tables for the eclipse black oil formulation\n using the oil component as brine and the gas component as co_2." << std::endl;
135 if (std::fabs(temperature-400.) > 100.0) {
136 std::cout <<
"T=" << temperature <<
" is outside the allowed range [300,500] Kelvin" << std::endl;
140 temperature_ = temperature;
147 std::ofstream black(
"blackoil_pvt");
149 black << std::fixed << std::showpoint;
151 const double pMin=150.e5;
152 const double pMax=400.e5;
153 const unsigned int np=11;
154 std::vector<double> pValue(np+1,0.0);
155 std::vector<double> rs(np+1,0.0);
157 pValue[0] = 101330.0;
162 for (
unsigned int i=0; i<np; ++i) {
163 pValue[i+1] = pMin + i*(pMax-pMin)/(np-1);
164 rs[i+1] =
R(pValue[i+1], z,
Liquid);
167 const unsigned int np_fill_in = 10;
168 const double dr = (rs[1] - rs[0])/(np_fill_in+1);
169 const double dp = (pValue[1] - pValue[0])/(np_fill_in+1);
170 rs.insert(rs.begin()+1,np_fill_in,0.0);
171 pValue.insert(pValue.begin()+1,np_fill_in,0.0);
172 for (
unsigned int i=1; i<=np_fill_in; ++i) {
173 rs[i] = rs[i-1] + dr;
174 pValue[i] = pValue[i-1] + dp;
179 black <<
"-- Rs Pbub Bo Vo\n";
180 black <<
"-- sm3/sm3 barsa rm3/sm3 cP\n";
183 for (
unsigned int i=0; i<np+np_fill_in; ++i) {
185 black << std::setw(14) << rs[i];
186 for (
unsigned int j=i; j<np+1+np_fill_in; ++j) {
188 if (j==i) black << std::setw(14) << pValue[j]*1.e-5 << std::setw(14) << 1.0-j*0.001 << std::setw(14) << 1.06499;
191 if (j>i) black << std::endl << std::setw(14) <<
' ';
192 black << std::setw(14) << pValue[j]*1.e-5
193 << std::setw(14) <<
B(pValue[j], z,
Liquid)
196 black <<
" /" << std::endl;
198 black <<
"/ " << std::endl;
204 black <<
"-- Pg Bg Vg\n";
205 black <<
"-- barsa rm3/sm3 cP\n";
207 for (
unsigned int i=0; i<np; ++i) {
210 black << std::setw(14) << pValue[i+np_fill_in+1]*1.e-5
211 << std::setw(14) <<
B(pValue[i+np_fill_in+1], z,
Vapour)
215 black <<
"/ " << std::endl;
219 black <<
"-- Pg Rv Bg Vg\n";
220 black <<
"-- barsa sm3/sm3 rm3/sm3 cP\n";
221 for (
unsigned int i=0; i<np; ++i) {
224 black << std::setw(14) << pValue[i+np_fill_in+1]*1.e-5
225 << std::setw(14) <<
R(pValue[i+np_fill_in+1], z,
Vapour)
226 << std::setw(14) <<
B(pValue[i+np_fill_in+1], z,
Vapour)
230 black << std::setw(14) <<
' '
231 << std::setw(14) <<
R(pValue[i+np_fill_in+1], z,
Vapour)
232 << std::setw(14) <<
B(pValue[i+np_fill_in+1], z,
Vapour)
234 <<
" /" << std::endl;
236 black <<
"/ " << std::endl;
238 std::cout <<
" Pvt tables for temperature=" << temperature <<
" Kelvin is written to file blackoil_pvt. " << std::endl;
239 std::cout <<
" NOTE that the file contains tables for both PVDG (dry gas) and PVTG (wet gas)." << std::endl;
244 return surfaceDensities_;
250 computeState(ss, surfvol[
Oil], surfvol[
Gas], press);
252 case Aqua:
return 1.0e-10;
253 case Liquid:
return ss.phaseViscosity[wPhase];
254 case Vapour:
return ss.phaseViscosity[nPhase];
263 computeState(ss, surfvol[
Oil], surfvol[
Gas], press);
265 case Aqua:
return 0.0;
266 case Liquid:
return ss.saturation;
267 case Vapour:
return 1.0 - ss.saturation;
276 computeState(ss, surfvol[
Oil], surfvol[
Gas], press);
278 case Aqua:
return 1.0;
279 case Liquid:
return surfaceDensities_[
Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10);
280 case Vapour:
return surfaceDensities_[
Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10);
287 const double dp = 100.;
288 return (
B(press+dp,surfvol,phase)-
B(press,surfvol,phase))/dp;
294 computeState(ss, surfvol[
Oil], surfvol[
Gas], press);
296 case Aqua:
return 0.0;
297 case Liquid:
return (ss.massfrac[wPhase][nComp]*surfaceDensities_[
Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[
Gas]+1.0e-10);
298 case Vapour:
return (ss.massfrac[nPhase][wComp]*surfaceDensities_[
Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[
Oil]+1.0e-10);
305 const double dp = 100.;
306 return (
R(press+dp,surfvol,phase)-
R(press,surfvol,phase))/dp;
310 const std::vector<CompVec>& surfvol,
311 std::vector<PhaseVec>& output)
const
313 int num = pressures.size();
316 for (
int i = 0; i < num; ++i) {
317 computeState(ss, surfvol[i][
Oil], surfvol[i][
Gas], pressures[i][
Liquid]);
318 output[i][
Aqua] = 1.0e-10;
319 output[i][
Liquid] = ss.phaseViscosity[wPhase];
320 output[i][
Vapour] = ss.phaseViscosity[nPhase];
325 const std::vector<CompVec>& surfvol,
326 std::vector<PhaseVec>& output)
const
328 int num = pressures.size();
331 for (
int i = 0; i < num; ++i) {
332 computeState(ss, surfvol[i][
Oil], surfvol[i][
Gas], pressures[i][
Liquid]);
333 output[i][
Aqua] = 1.0;
334 output[i][
Liquid] = surfaceDensities_[
Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10);
335 output[i][
Vapour] = surfaceDensities_[
Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10);
340 const std::vector<CompVec>& surfvol,
341 std::vector<PhaseVec>& output_B,
342 std::vector<PhaseVec>& output_dBdp)
const
344 int num = pressures.size();
345 output_B.resize(num);
346 output_dBdp.resize(num);
348 const double dp = 100.;
349 for (
int i = 0; i < num; ++i) {
350 computeState(ss, surfvol[i][
Oil], surfvol[i][
Gas], pressures[i][
Liquid]);
351 output_B[i][
Aqua] = 1.0;
352 output_B[i][
Liquid] = surfaceDensities_[
Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10);
353 output_B[i][
Vapour] = surfaceDensities_[
Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10);
354 computeState(ss, surfvol[i][
Oil], surfvol[i][
Gas], pressures[i][
Liquid]+dp);
355 output_dBdp[i][
Aqua] = 0.0;
356 output_dBdp[i][
Liquid] = (surfaceDensities_[
Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10) - output_B[i][
Liquid])/dp;
357 output_dBdp[i][
Vapour] = (surfaceDensities_[
Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10) - output_B[i][
Vapour])/dp;
362 const std::vector<CompVec>& surfvol,
363 std::vector<PhaseVec>& output)
const
365 int num = pressures.size();
368 for (
int i = 0; i < num; ++i) {
369 computeState(ss, surfvol[i][
Oil], surfvol[i][
Gas], pressures[i][
Liquid]);
370 output[i][
Aqua] = 0.0;
371 output[i][
Liquid] = (ss.massfrac[wPhase][nComp]*surfaceDensities_[
Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[
Gas]+1.0e-10);
372 output[i][
Vapour] = (ss.massfrac[nPhase][wComp]*surfaceDensities_[
Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[
Oil]+1.0e-10);
377 const std::vector<CompVec>& surfvol,
378 std::vector<PhaseVec>& output_R,
379 std::vector<PhaseVec>& output_dRdp)
const
381 int num = pressures.size();
382 output_R.resize(num);
383 output_dRdp.resize(num);
385 const double dp = 100.;
386 for (
int i = 0; i < num; ++i) {
387 computeState(ss, surfvol[i][
Oil], surfvol[i][
Gas], pressures[i][
Liquid]);
388 output_R[i][
Aqua] = 0.0;
389 output_R[i][
Liquid] = (ss.massfrac[wPhase][nComp]*surfaceDensities_[
Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[
Gas]+1.0e-10);
390 output_R[i][
Vapour] = (ss.massfrac[nPhase][wComp]*surfaceDensities_[
Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[
Oil]+1.0e-10);
391 computeState(ss, surfvol[i][
Oil], surfvol[i][
Gas], pressures[i][
Liquid]+dp);
392 output_dRdp[i][
Aqua] = 0.0;
393 output_dRdp[i][
Liquid] = ((ss.massfrac[wPhase][nComp]*surfaceDensities_[
Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[
Gas]+1.0e-10) - output_R[i][
Liquid])/dp;
394 output_dRdp[i][
Vapour] = ((ss.massfrac[nPhase][wComp]*surfaceDensities_[
Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[
Oil]+1.0e-10) - output_R[i][
Vapour])/dp;
398void BlackoilCo2PVT::computeState(BlackoilCo2PVT::SubState& ss,
double zBrine,
double zCO2,
double pressure)
const
402 fluidState.setTemperature(temperature_);
403 fluidState.setPressure(wPhase, pressure);
404 fluidState.setPressure(nPhase, pressure);
406 double massH20 = surfaceDensities_[
Oil]*zBrine;
407 double massCO2 = surfaceDensities_[
Gas]*zCO2;
410 typename FluidSystem::ParameterCache paramCache;
411 typedef Opm::MiscibleMultiPhaseComposition<double,
FluidSystem> MMPC;
412 MMPC::solve(fluidState, paramCache,
false,
false);
413 ss.density[wPhase] = fluidState.density(wPhase);
414 ss.density[nPhase] = fluidState.density(nPhase);
415 ss.massfrac[wPhase][nComp] = fluidState.massFraction(wPhase, nComp);
416 ss.massfrac[nPhase][wComp] = fluidState.massFraction(nPhase, wComp);
417 ss.massfrac[wPhase][wComp] = 1.0 - ss.massfrac[wPhase][nComp];
418 ss.massfrac[nPhase][nComp] = 1.0 - ss.massfrac[nPhase][wComp];
420 double detX = ss.massfrac[wPhase][wComp]*ss.massfrac[nPhase][nComp]-ss.massfrac[wPhase][nComp]*ss.massfrac[nPhase][wComp];
421 ss.phaseVolume[wPhase] = (massH20*ss.massfrac[nPhase][nComp] - massCO2*ss.massfrac[nPhase][wComp])/(ss.density[wPhase]*detX);
422 ss.phaseVolume[nPhase] = (massCO2*ss.massfrac[wPhase][wComp] - massH20*ss.massfrac[wPhase][nComp])/(ss.density[nPhase]*detX);
425 if (ss.phaseVolume[wPhase] > 0.0 && ss.phaseVolume[nPhase] > 0.0) {
426 ss.saturation = ss.phaseVolume[wPhase]/(ss.phaseVolume[wPhase]+ss.phaseVolume[nPhase]);
427 fluidState.setSaturation(wPhase, ss.saturation);
428 fluidState.setSaturation(nPhase, 1.0 - ss.saturation);
430 else if (ss.phaseVolume[wPhase] <= 0.0) {
433 ss.massfrac[nPhase][nComp] = massCO2/(massCO2+massH20);
434 ss.massfrac[nPhase][wComp] = 1.0 - ss.massfrac[nPhase][nComp];
435 double M1 = FluidSystem::molarMass(wComp);
436 double M2 = FluidSystem::molarMass(nComp);
437 double avgMolarMass = M1*M2/(M2 + ss.massfrac[nPhase][nComp]*(M1 - M2));
438 fluidState.setMoleFraction(nPhase, nComp, ss.massfrac[nPhase][nComp]*avgMolarMass/M2);
439 fluidState.setMoleFraction(nPhase, wComp, ss.massfrac[nPhase][wComp]*avgMolarMass/M1);
440 ss.density[nPhase] = brineCo2_.density(fluidState, paramCache, nPhase);
441 fluidState.setDensity(nPhase, ss.density[nPhase]);
442 ss.phaseVolume[nPhase] = (massH20+massCO2)/ss.density[nPhase];
443 fluidState.setSaturation(nPhase, 1.0 - ss.saturation);
445 paramCache.updatePhase(fluidState, nPhase);
446 typedef Opm::ComputeFromReferencePhase<double,
FluidSystem> CFRP;
447 CFRP::solve(fluidState,
452 ss.massfrac[wPhase][wComp] = fluidState.massFraction(wPhase, wComp);
453 ss.massfrac[wPhase][nComp] = fluidState.massFraction(wPhase, nComp);
454 ss.density[wPhase] = fluidState.density(wPhase);
455 ss.phaseVolume[wPhase] = 0.0;
456 fluidState.setSaturation(wPhase, ss.saturation);
458 else if (ss.phaseVolume[nPhase] <= 0.0) {
461 ss.massfrac[wPhase][wComp] = massH20/(massCO2+massH20);
462 ss.massfrac[wPhase][nComp] = 1.0 - ss.massfrac[wPhase][wComp];
463 double M1 = FluidSystem::molarMass(wComp);
464 double M2 = FluidSystem::molarMass(nComp);
465 double avgMolarMass = M1*M2/(M2 + ss.massfrac[wPhase][nComp]*(M1 - M2));
466 fluidState.setMoleFraction(wPhase, nComp, ss.massfrac[wPhase][nComp]*avgMolarMass/M2);
467 fluidState.setMoleFraction(wPhase, wComp, ss.massfrac[wPhase][wComp]*avgMolarMass/M1);
468 ss.density[wPhase] = brineCo2_.density(fluidState, paramCache, wPhase);
469 fluidState.setDensity(wPhase, ss.density[wPhase]);
470 ss.phaseVolume[wPhase] = (massH20+massCO2)/ss.density[wPhase];
471 fluidState.setSaturation(wPhase, ss.saturation);
473 paramCache.updatePhase(fluidState, nPhase);
474 typedef ComputeFromReferencePhase<double,
FluidSystem> CFRP;
475 CFRP::solve(fluidState,
480 ss.massfrac[nPhase][nComp] = fluidState.massFraction(nPhase, nComp);
481 ss.massfrac[nPhase][wComp] = fluidState.massFraction(nPhase, wComp);
482 ss.density[nPhase] = fluidState.density(nPhase);
483 ss.phaseVolume[nPhase] = 0.0;
484 fluidState.setSaturation(nPhase, 1.0 - ss.saturation);
487 ss.phaseViscosity[wPhase] = brineCo2_.viscosity(fluidState, paramCache, wPhase);
488 ss.phaseViscosity[nPhase] = brineCo2_.viscosity(fluidState, paramCache, nPhase);
Provides the class with the tabulated values of CO2 for the benchmark3 problem.
Definition: BlackoilCo2PVT.hpp:48
double getSaturation(double press, const CompVec &surfvol, PhaseIndex phase) const
Definition: BlackoilCo2PVT.hpp:260
double dBdp(double press, const CompVec &surfvol, PhaseIndex phase) const
Definition: BlackoilCo2PVT.hpp:285
CompVec surfaceDensities() const
Definition: BlackoilCo2PVT.hpp:242
double dRdp(double press, const CompVec &surfvol, PhaseIndex phase) const
Definition: BlackoilCo2PVT.hpp:303
Opm::CompositionalFluidState< double, FluidSystem > CompositionalFluidState
Definition: BlackoilCo2PVT.hpp:52
double B(double press, const CompVec &surfvol, PhaseIndex phase) const
Definition: BlackoilCo2PVT.hpp:273
void generateBlackOilTables(double temperature)
Definition: BlackoilCo2PVT.hpp:132
void init(Opm::DeckConstPtr deck)
Definition: BlackoilCo2PVT.hpp:121
double R(double press, const CompVec &surfvol, PhaseIndex phase) const
Definition: BlackoilCo2PVT.hpp:291
double getViscosity(double press, const CompVec &surfvol, PhaseIndex phase) const
Definition: BlackoilCo2PVT.hpp:247
Opm::FluidSystems::BrineCO2< double, Opm::Benchmark3::CO2Tables > FluidSystem
Definition: BlackoilCo2PVT.hpp:51
Definition: BlackoilDefs.hpp:31
Dune::FieldVector< Scalar, numComponents > CompVec
Definition: BlackoilDefs.hpp:40
PhaseIndex
Definition: BlackoilDefs.hpp:37
@ Aqua
Definition: BlackoilDefs.hpp:37
@ Vapour
Definition: BlackoilDefs.hpp:37
@ Liquid
Definition: BlackoilDefs.hpp:37
@ Gas
Definition: BlackoilDefs.hpp:36
@ Oil
Definition: BlackoilDefs.hpp:36
@ Water
Definition: BlackoilDefs.hpp:36
Definition: BlackoilFluid.hpp:32
BlackoilCo2PVT BlackoilPVT
Definition: BlackoilCo2PVT.hpp:494