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>
35 #include <opm/porsol/blackoil/fluid/MiscibilityProps.hpp>
36 #include <opm/porsol/blackoil/fluid/BlackoilDefs.hpp>
38 #include <opm/parser/eclipse/Deck/Deck.hpp>
53 typedef Opm::CompositionalFluidState<double, FluidSystem> CompositionalFluidState;
55 void init(
const Opm::Deck& deck);
57 void generateBlackOilTables(
double temperature);
59 double getViscosity(
double press,
60 const CompVec& surfvol,
61 PhaseIndex phase)
const;
62 CompVec surfaceDensities()
const;
63 double getSaturation(
double press,
64 const CompVec& surfvol,
65 PhaseIndex phase)
const;
66 double B (
double press,
67 const CompVec& surfvol,
68 PhaseIndex phase)
const;
69 double dBdp(
double press,
70 const CompVec& surfvol,
71 PhaseIndex phase)
const;
72 double R (
double press,
73 const CompVec& surfvol,
74 PhaseIndex phase)
const;
75 double dRdp(
double press,
76 const CompVec& surfvol,
77 PhaseIndex phase)
const;
79 void getViscosity(
const std::vector<PhaseVec>& pressures,
80 const std::vector<CompVec>& surfvol,
81 std::vector<PhaseVec>& output)
const;
82 void B(
const std::vector<PhaseVec>& pressures,
83 const std::vector<CompVec>& surfvol,
84 std::vector<PhaseVec>& output)
const;
85 void dBdp(
const std::vector<PhaseVec>& pressures,
86 const std::vector<CompVec>& surfvol,
87 std::vector<PhaseVec>& output_B,
88 std::vector<PhaseVec>& output_dBdp)
const;
89 void R(
const std::vector<PhaseVec>& pressures,
90 const std::vector<CompVec>& surfvol,
91 std::vector<PhaseVec>& output)
const;
92 void dRdp(
const std::vector<PhaseVec>& pressures,
93 const std::vector<CompVec>& surfvol,
94 std::vector<PhaseVec>& output_R,
95 std::vector<PhaseVec>& output_dRdp)
const;
98 CompVec surfaceDensities_;
99 FluidSystem brineCo2_;
103 Dune::FieldVector<double,2> density;
104 Dune::FieldMatrix<double,2,2> massfrac;
105 Dune::FieldVector<double,2> phaseVolume;
106 Dune::FieldVector<double,2> phaseViscosity;
109 void computeState(SubState& ss,
double zBrine,
double zCO2,
double pressure)
const;
111 wPhase = FluidSystem::liquidPhaseIdx,
112 nPhase = FluidSystem::gasPhaseIdx,
114 wComp = FluidSystem::BrineIdx,
115 nComp = FluidSystem::CO2Idx
122 void BlackoilCo2PVT::init(
const Opm::Deck& )
124 surfaceDensities_[Water] = 1000.;
125 surfaceDensities_[Gas] = 2.0;
126 surfaceDensities_[Oil] = 1000.;
133 void BlackoilCo2PVT::generateBlackOilTables(
double temperature)
135 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;
136 if (std::fabs(temperature-400.) > 100.0) {
137 std::cout <<
"T=" << temperature <<
" is outside the allowed range [300,500] Kelvin" << std::endl;
141 temperature_ = temperature;
148 std::ofstream black(
"blackoil_pvt");
150 black << std::fixed << std::showpoint;
152 const double pMin=150.e5;
153 const double pMax=400.e5;
154 const unsigned int np=11;
155 std::vector<double> pValue(np+1,0.0);
156 std::vector<double> rs(np+1,0.0);
158 pValue[0] = 101330.0;
163 for (
unsigned int i=0; i<np; ++i) {
164 pValue[i+1] = pMin + i*(pMax-pMin)/(np-1);
165 rs[i+1] = R(pValue[i+1], z, Liquid);
168 const unsigned int np_fill_in = 10;
169 const double dr = (rs[1] - rs[0])/(np_fill_in+1);
170 const double dp = (pValue[1] - pValue[0])/(np_fill_in+1);
171 rs.insert(rs.begin()+1,np_fill_in,0.0);
172 pValue.insert(pValue.begin()+1,np_fill_in,0.0);
173 for (
unsigned int i=1; i<=np_fill_in; ++i) {
174 rs[i] = rs[i-1] + dr;
175 pValue[i] = pValue[i-1] + dp;
180 black <<
"-- Rs Pbub Bo Vo\n";
181 black <<
"-- sm3/sm3 barsa rm3/sm3 cP\n";
184 for (
unsigned int i=0; i<np+np_fill_in; ++i) {
186 black << std::setw(14) << rs[i];
187 for (
unsigned int j=i; j<np+1+np_fill_in; ++j) {
189 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;
192 if (j>i) black << std::endl << std::setw(14) <<
' ';
193 black << std::setw(14) << pValue[j]*1.e-5
194 << std::setw(14) << B(pValue[j], z, Liquid)
195 << std::setw(14) << getViscosity(pValue[j], z, Liquid)*1.e3;
197 black <<
" /" << std::endl;
199 black <<
"/ " << std::endl;
205 black <<
"-- Pg Bg Vg\n";
206 black <<
"-- barsa rm3/sm3 cP\n";
208 for (
unsigned int i=0; i<np; ++i) {
211 black << std::setw(14) << pValue[i+np_fill_in+1]*1.e-5
212 << std::setw(14) << B(pValue[i+np_fill_in+1], z, Vapour)
213 << std::setw(14) << getViscosity(pValue[i+np_fill_in+1], z, Vapour)*1.e3
216 black <<
"/ " << std::endl;
220 black <<
"-- Pg Rv Bg Vg\n";
221 black <<
"-- barsa sm3/sm3 rm3/sm3 cP\n";
222 for (
unsigned int i=0; i<np; ++i) {
225 black << std::setw(14) << pValue[i+np_fill_in+1]*1.e-5
226 << std::setw(14) << R(pValue[i+np_fill_in+1], z, Vapour)
227 << std::setw(14) << B(pValue[i+np_fill_in+1], z, Vapour)
228 << std::setw(14) << getViscosity(pValue[i+np_fill_in+1], z, Vapour)*1.e3
231 black << std::setw(14) <<
' '
232 << std::setw(14) << R(pValue[i+np_fill_in+1], z, Vapour)
233 << std::setw(14) << B(pValue[i+np_fill_in+1], z, Vapour)
234 << std::setw(14) << getViscosity(pValue[i+np_fill_in+1], z, Vapour)*1.e3
235 <<
" /" << std::endl;
237 black <<
"/ " << std::endl;
239 std::cout <<
" Pvt tables for temperature=" << temperature <<
" Kelvin is written to file blackoil_pvt. " << std::endl;
240 std::cout <<
" NOTE that the file contains tables for both PVDG (dry gas) and PVTG (wet gas)." << std::endl;
243 BlackoilCo2PVT::CompVec BlackoilCo2PVT::surfaceDensities()
const
245 return surfaceDensities_;
248 double BlackoilCo2PVT::getViscosity(
double press,
const CompVec& surfvol, PhaseIndex phase)
const
251 computeState(ss, surfvol[Oil], surfvol[Gas], press);
253 case Aqua:
return 1.0e-10;
254 case Liquid:
return ss.phaseViscosity[wPhase];
255 case Vapour:
return ss.phaseViscosity[nPhase];
261 double BlackoilCo2PVT::getSaturation(
double press,
const CompVec& surfvol, PhaseIndex phase)
const
264 computeState(ss, surfvol[Oil], surfvol[Gas], press);
266 case Aqua:
return 0.0;
267 case Liquid:
return ss.saturation;
268 case Vapour:
return 1.0 - ss.saturation;
274 double BlackoilCo2PVT::B(
double press,
const CompVec& surfvol, PhaseIndex phase)
const
277 computeState(ss, surfvol[Oil], surfvol[Gas], press);
279 case Aqua:
return 1.0;
280 case Liquid:
return surfaceDensities_[Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10);
281 case Vapour:
return surfaceDensities_[Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10);
286 double BlackoilCo2PVT::dBdp(
double press,
const CompVec& surfvol, PhaseIndex phase)
const
288 const double dp = 100.;
289 return (B(press+dp,surfvol,phase)-B(press,surfvol,phase))/dp;
292 double BlackoilCo2PVT::R(
double press,
const CompVec& surfvol, PhaseIndex phase)
const
295 computeState(ss, surfvol[Oil], surfvol[Gas], press);
297 case Aqua:
return 0.0;
298 case Liquid:
return (ss.massfrac[wPhase][nComp]*surfaceDensities_[Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[Gas]+1.0e-10);
299 case Vapour:
return (ss.massfrac[nPhase][wComp]*surfaceDensities_[Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[Oil]+1.0e-10);
304 double BlackoilCo2PVT::dRdp(
double press,
const CompVec& surfvol, PhaseIndex phase)
const
306 const double dp = 100.;
307 return (R(press+dp,surfvol,phase)-R(press,surfvol,phase))/dp;
310 void BlackoilCo2PVT::getViscosity(
const std::vector<PhaseVec>& pressures,
311 const std::vector<CompVec>& surfvol,
312 std::vector<PhaseVec>& output)
const
314 int num = pressures.size();
317 for (
int i = 0; i < num; ++i) {
318 computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]);
319 output[i][Aqua] = 1.0e-10;
320 output[i][Liquid] = ss.phaseViscosity[wPhase];
321 output[i][Vapour] = ss.phaseViscosity[nPhase];
325 void BlackoilCo2PVT::B(
const std::vector<PhaseVec>& pressures,
326 const std::vector<CompVec>& surfvol,
327 std::vector<PhaseVec>& output)
const
329 int num = pressures.size();
332 for (
int i = 0; i < num; ++i) {
333 computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]);
334 output[i][Aqua] = 1.0;
335 output[i][Liquid] = surfaceDensities_[Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10);
336 output[i][Vapour] = surfaceDensities_[Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10);
340 void BlackoilCo2PVT::dBdp(
const std::vector<PhaseVec>& pressures,
341 const std::vector<CompVec>& surfvol,
342 std::vector<PhaseVec>& output_B,
343 std::vector<PhaseVec>& output_dBdp)
const
345 int num = pressures.size();
346 output_B.resize(num);
347 output_dBdp.resize(num);
349 const double dp = 100.;
350 for (
int i = 0; i < num; ++i) {
351 computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]);
352 output_B[i][Aqua] = 1.0;
353 output_B[i][Liquid] = surfaceDensities_[Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10);
354 output_B[i][Vapour] = surfaceDensities_[Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10);
355 computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]+dp);
356 output_dBdp[i][Aqua] = 0.0;
357 output_dBdp[i][Liquid] = (surfaceDensities_[Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10) - output_B[i][Liquid])/dp;
358 output_dBdp[i][Vapour] = (surfaceDensities_[Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10) - output_B[i][Vapour])/dp;
362 void BlackoilCo2PVT::R(
const std::vector<PhaseVec>& pressures,
363 const std::vector<CompVec>& surfvol,
364 std::vector<PhaseVec>& output)
const
366 int num = pressures.size();
369 for (
int i = 0; i < num; ++i) {
370 computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]);
371 output[i][Aqua] = 0.0;
372 output[i][Liquid] = (ss.massfrac[wPhase][nComp]*surfaceDensities_[Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[Gas]+1.0e-10);
373 output[i][Vapour] = (ss.massfrac[nPhase][wComp]*surfaceDensities_[Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[Oil]+1.0e-10);
377 void BlackoilCo2PVT::dRdp(
const std::vector<PhaseVec>& pressures,
378 const std::vector<CompVec>& surfvol,
379 std::vector<PhaseVec>& output_R,
380 std::vector<PhaseVec>& output_dRdp)
const
382 int num = pressures.size();
383 output_R.resize(num);
384 output_dRdp.resize(num);
386 const double dp = 100.;
387 for (
int i = 0; i < num; ++i) {
388 computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]);
389 output_R[i][Aqua] = 0.0;
390 output_R[i][Liquid] = (ss.massfrac[wPhase][nComp]*surfaceDensities_[Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[Gas]+1.0e-10);
391 output_R[i][Vapour] = (ss.massfrac[nPhase][wComp]*surfaceDensities_[Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[Oil]+1.0e-10);
392 computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]+dp);
393 output_dRdp[i][Aqua] = 0.0;
394 output_dRdp[i][Liquid] = ((ss.massfrac[wPhase][nComp]*surfaceDensities_[Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[Gas]+1.0e-10) - output_R[i][Liquid])/dp;
395 output_dRdp[i][Vapour] = ((ss.massfrac[nPhase][wComp]*surfaceDensities_[Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[Oil]+1.0e-10) - output_R[i][Vapour])/dp;
399 void BlackoilCo2PVT::computeState(BlackoilCo2PVT::SubState& ss,
double zBrine,
double zCO2,
double pressure)
const
402 CompositionalFluidState fluidState;
403 fluidState.setTemperature(temperature_);
404 fluidState.setPressure(wPhase, pressure);
405 fluidState.setPressure(nPhase, pressure);
407 double massH20 = surfaceDensities_[Oil]*zBrine;
408 double massCO2 = surfaceDensities_[Gas]*zCO2;
411 FluidSystem::ParameterCache<double> paramCache;
412 typedef Opm::MiscibleMultiPhaseComposition<double, FluidSystem> MMPC;
413 MMPC::solve(fluidState, paramCache,
false,
false);
414 ss.density[wPhase] = fluidState.density(wPhase);
415 ss.density[nPhase] = fluidState.density(nPhase);
416 ss.massfrac[wPhase][nComp] = fluidState.massFraction(wPhase, nComp);
417 ss.massfrac[nPhase][wComp] = fluidState.massFraction(nPhase, wComp);
418 ss.massfrac[wPhase][wComp] = 1.0 - ss.massfrac[wPhase][nComp];
419 ss.massfrac[nPhase][nComp] = 1.0 - ss.massfrac[nPhase][wComp];
421 double detX = ss.massfrac[wPhase][wComp]*ss.massfrac[nPhase][nComp]-ss.massfrac[wPhase][nComp]*ss.massfrac[nPhase][wComp];
422 ss.phaseVolume[wPhase] = (massH20*ss.massfrac[nPhase][nComp] - massCO2*ss.massfrac[nPhase][wComp])/(ss.density[wPhase]*detX);
423 ss.phaseVolume[nPhase] = (massCO2*ss.massfrac[wPhase][wComp] - massH20*ss.massfrac[wPhase][nComp])/(ss.density[nPhase]*detX);
426 if (ss.phaseVolume[wPhase] > 0.0 && ss.phaseVolume[nPhase] > 0.0) {
427 ss.saturation = ss.phaseVolume[wPhase]/(ss.phaseVolume[wPhase]+ss.phaseVolume[nPhase]);
428 fluidState.setSaturation(wPhase, ss.saturation);
429 fluidState.setSaturation(nPhase, 1.0 - ss.saturation);
431 else if (ss.phaseVolume[wPhase] <= 0.0) {
434 ss.massfrac[nPhase][nComp] = massCO2/(massCO2+massH20);
435 ss.massfrac[nPhase][wComp] = 1.0 - ss.massfrac[nPhase][nComp];
436 double M1 = FluidSystem::molarMass(wComp);
437 double M2 = FluidSystem::molarMass(nComp);
438 double avgMolarMass = M1*M2/(M2 + ss.massfrac[nPhase][nComp]*(M1 - M2));
439 fluidState.setMoleFraction(nPhase, nComp, ss.massfrac[nPhase][nComp]*avgMolarMass/M2);
440 fluidState.setMoleFraction(nPhase, wComp, ss.massfrac[nPhase][wComp]*avgMolarMass/M1);
441 ss.density[nPhase] = brineCo2_.density(fluidState, paramCache, nPhase);
442 fluidState.setDensity(nPhase, ss.density[nPhase]);
443 ss.phaseVolume[nPhase] = (massH20+massCO2)/ss.density[nPhase];
444 fluidState.setSaturation(nPhase, 1.0 - ss.saturation);
446 paramCache.updatePhase(fluidState, nPhase);
447 typedef Opm::ComputeFromReferencePhase<double, FluidSystem> CFRP;
448 CFRP::solve(fluidState,
453 ss.massfrac[wPhase][wComp] = fluidState.massFraction(wPhase, wComp);
454 ss.massfrac[wPhase][nComp] = fluidState.massFraction(wPhase, nComp);
455 ss.density[wPhase] = fluidState.density(wPhase);
456 ss.phaseVolume[wPhase] = 0.0;
457 fluidState.setSaturation(wPhase, ss.saturation);
459 else if (ss.phaseVolume[nPhase] <= 0.0) {
462 ss.massfrac[wPhase][wComp] = massH20/(massCO2+massH20);
463 ss.massfrac[wPhase][nComp] = 1.0 - ss.massfrac[wPhase][wComp];
464 double M1 = FluidSystem::molarMass(wComp);
465 double M2 = FluidSystem::molarMass(nComp);
466 double avgMolarMass = M1*M2/(M2 + ss.massfrac[wPhase][nComp]*(M1 - M2));
467 fluidState.setMoleFraction(wPhase, nComp, ss.massfrac[wPhase][nComp]*avgMolarMass/M2);
468 fluidState.setMoleFraction(wPhase, wComp, ss.massfrac[wPhase][wComp]*avgMolarMass/M1);
469 ss.density[wPhase] = brineCo2_.density(fluidState, paramCache, wPhase);
470 fluidState.setDensity(wPhase, ss.density[wPhase]);
471 ss.phaseVolume[wPhase] = (massH20+massCO2)/ss.density[wPhase];
472 fluidState.setSaturation(wPhase, ss.saturation);
474 paramCache.updatePhase(fluidState, nPhase);
475 typedef ComputeFromReferencePhase<double, FluidSystem> CFRP;
476 CFRP::solve(fluidState,
481 ss.massfrac[nPhase][nComp] = fluidState.massFraction(nPhase, nComp);
482 ss.massfrac[nPhase][wComp] = fluidState.massFraction(nPhase, wComp);
483 ss.density[nPhase] = fluidState.density(nPhase);
484 ss.phaseVolume[nPhase] = 0.0;
485 fluidState.setSaturation(nPhase, 1.0 - ss.saturation);
488 ss.phaseViscosity[wPhase] = brineCo2_.viscosity(fluidState, paramCache, wPhase);
489 ss.phaseViscosity[nPhase] = brineCo2_.viscosity(fluidState, paramCache, nPhase);
495 typedef BlackoilCo2PVT BlackoilPVT;
500 #endif // OPM_BLACKOILCO2PVT_HEADER_INCLUDED
Provides the class with the tabulated values of CO2 for the benchmark3 problem.
Definition: BlackoilDefs.hpp:33
Definition: benchmark3co2tables.hh:25261
Definition: BlackoilCo2PVT.hpp:48