Added pvt functionality for wetgas

The pvt interface is extended to handle wet-gas systems:
1. rvSat is added as a function in the PVT interface
2. SinglePvtLiveGas computes the pvt values and its derivatives
3. The old rbub variable is changed to rsSat for clearity
4. The new interface is tested in test_blackoilfluid with data from
liveoil.DATA and wetgas.DATA
This commit is contained in:
Tor Harald Sandve committed 2014-01-10 16:07:02 +01:00
1 parent 15211bc445
commit 598c7a0cb8
16 files changed
+604 -165

No files matched your search

+36
View File
@@ -0,0 +1,36 @@
OIL
GAS
WATER
FIELD
PVTO
-- Rs Pbub Bo Vo
.0 14.7 1.0000 1.20 /
.165 400. 1.0120 1.17 /
.335 800. 1.0255 1.14 /
.500 1200. 1.0380 1.11 /
.665 1600. 1.0510 1.08 /
.828 2000. 1.0630 1.06 /
.985 2400. 1.0750 1.03 /
1.130 2800. 1.0870 1.00 /
1.270 3200. 1.0985 .98 /
1.390 3600. 1.1100 .95 /
1.500 4000. 1.1200 .94
5000. 1.1189 .94 /
/
PVDG
-- Pg Bg Vg
14.7 178.08 .0125
400. 5.4777 .0130
800. 2.7392 .0135
1200. 1.8198 .0140
1600. 1.3648 .0145
2000. 1.0957 .0150
2400. 0.9099 .0155
2800. 0.7799 .0160
3200. 0.6871 .0165
3600. 0.6035 .0170
4000. 0.5432 .0175 /
+224 -103
View File
@@ -28,24 +28,18 @@
using namespace Opm;
using namespace std;
BOOST_AUTO_TEST_CASE(test_blackoilfluid)
{
std::vector<std::shared_ptr<SinglePvtInterface> > getProps(const EclipseGridParser deck,PhaseUsage phase_usage_){
// read eclipse deck
const string filename = "testFluid.DATA";
cout << "Reading deck: " << filename << endl;
const EclipseGridParser deck (filename);
// setup pvt interface
std::vector<std::shared_ptr<SinglePvtInterface> > props_;
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
int samples = 0;
std::vector<std::shared_ptr<SinglePvtInterface> > props_;
// Set the properties.
props_.resize(phase_usage_.num_phases);
// Water PVT
if (phase_usage_.phase_used[Aqua]) {
if (deck.hasField("PVTW")) {
@@ -55,6 +49,7 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
}
}
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
if (deck.hasField("PVDO")) {
@@ -87,6 +82,146 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
}
}
return props_;
}
void testmu(const double reltol, int n, int np, std::vector<double> p, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<SinglePvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
{
std::vector<double> mu(n);
std::vector<double> dmudp(n);
std::vector<double> dmudr(n);
std::vector<double> mu_new(n);
double dmudp_diff;
double dmudr_diff;
double dmudp_diff_u;
double dmudr_diff_u;
// test mu
for (int phase = 0; phase < np; ++phase) {
props_[phase]->mu(n, &p[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
props_[phase]->mu(n, &p[0], &z[0], &mu[0]);
dmudp_diff = (mu_new[1]-mu_new[0])/(p[1]-p[0]);
dmudr_diff = (mu_new[2]-mu_new[0])/(r[2]-r[0]);
dmudp_diff_u = (mu_new[4]-mu_new[3])/(p[4]-p[3]);
dmudr_diff_u = (mu_new[5]-mu_new[3])/(r[5]-r[3]);
for (int i = 0; i < n; ++i){
BOOST_CHECK_CLOSE(mu_new[i],mu[i],reltol);
}
// saturated case
BOOST_CHECK_CLOSE(dmudp_diff,dmudp[0],reltol);
BOOST_CHECK_CLOSE(dmudr_diff,dmudr[0],reltol);
// unsaturated case
BOOST_CHECK_CLOSE(dmudp_diff_u,dmudp[3],reltol);
BOOST_CHECK_CLOSE(dmudr_diff_u,dmudr[3],reltol);
}
}
void testb(const double reltol, int n, int np, std::vector<double> p, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<SinglePvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
{
// test b
std::vector<double> b(n);
std::vector<double> B(n);
std::vector<double> invB(n);
std::vector<double> dinvBdp(n);
std::vector<double> dBdp(n);
std::vector<double> dbdr(n);
std::vector<double> dbdp(n);
double dbdp_diff;
double dbdr_diff;
double dbdp_diff_u;
double dbdr_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase]->b(n, &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]);
props_[phase]->dBdp(n, &p[0], &z[0], &B[0], &dBdp[0]);
dbdp_diff = (b[1]-b[0])/(p[1]-p[0]);
dbdr_diff = (b[2]-b[0])/(r[2]-r[0]);
dbdp_diff_u = (b[4]-b[3])/(p[4]-p[3]);
dbdr_diff_u = (b[5]-b[3])/(r[5]-r[3]);
for (int i = 0; i < n; ++i){
invB[i] = 1/B[i];
dinvBdp[i] = -1/pow(B[i],2) * dBdp[i];
}
for (int i = 0; i < n; ++i){
BOOST_CHECK_CLOSE(invB[i],b[i] , reltol);
BOOST_CHECK_CLOSE(dinvBdp[i],dbdp[i] , reltol);
}
// saturated case
BOOST_CHECK_CLOSE(dbdp_diff,dbdp[0], reltol);
BOOST_CHECK_CLOSE(dbdr_diff,dbdr[0], reltol);
// unsaturated case
BOOST_CHECK_CLOSE(dbdp_diff_u,dbdp[3], reltol);
BOOST_CHECK_CLOSE(dbdr_diff_u,dbdr[3], reltol);
}
}
void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
// test bublepoint pressure
std::vector<double> rs(n);
std::vector<double> drsdp(n);
double drsdp_diff;
double drsdp_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase] ->rsSat(n, &p[0], &rs[0], &drsdp[0]);
drsdp_diff = (rs[1]-rs[0])/(p[1]-p[0]);
drsdp_diff_u = (rs[4]-rs[3])/(p[4]-p[3]);
// saturated case
BOOST_CHECK_CLOSE(drsdp_diff,drsdp[0], reltol);
// unsaturad case
BOOST_CHECK_CLOSE(drsdp_diff_u,drsdp[3], reltol);
}
}
void testrvSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
// test rv saturated
std::vector<double> rv(n);
std::vector<double> drvdp(n);
double drvdp_diff;
double drvdp_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase] ->rvSat(n, &p[0], &rv[0], &drvdp[0]);
drvdp_diff = (rv[1]-rv[0])/(p[1]-p[0]);
drvdp_diff_u = (rv[4]-rv[3])/(p[4]-p[3]);
// saturated case
BOOST_CHECK_CLOSE(drvdp_diff,drvdp[0], reltol);
// unsaturad case
BOOST_CHECK_CLOSE(drvdp_diff_u,drvdp[3], reltol);
}
}
BOOST_AUTO_TEST_CASE(test_liveoil)
{
// read eclipse deck
const string filename = "liveoil.DATA";
cout << "Reading deck: " << filename << endl;
const EclipseGridParser deck (filename);
// setup pvt interface
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(deck,phase_usage_);
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
// approximation of the derivatives.
const int n = 6;
@@ -97,25 +232,16 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
std::vector<double> p(n);
std::vector<double> r(n);
std::vector<Opm::PhasePresence> condition(n);
std::vector<PhasePresence> condition(n);
std::vector<double> z(n * np);
std::vector<double> mu(n);
std::vector<double> dmudp(n);
std::vector<double> dmudr(n);
std::vector<double> mu_new(n);
double dmudp_diff;
double dmudr_diff;
double dmudp_diff_u;
double dmudr_diff_u;
// Used for forward difference calculations
const double h_p = 100000;
const double h_p = 1e4;
const double h_r = 1;
// saturated
p[0] = 10000000;
p[0] = 1e7;
p[1] = p[0] + h_p;
p[2] = p[0];
@@ -144,90 +270,85 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
}
// test mu
for (int phase = 1; phase < 2; ++phase) {
props_[phase]->mu(n, &p[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
props_[phase]->mu(n, &p[0], &z[0], &mu[0]);
dmudp_diff = (mu_new[1]-mu_new[0])/h_p;
dmudr_diff = (mu_new[2]-mu_new[0])/h_r;
dmudp_diff_u = (mu_new[4]-mu_new[3])/h_p;
dmudr_diff_u = (mu_new[5]-mu_new[3])/h_r;
for (int i = 0; i < n; ++i){
BOOST_CHECK_CLOSE(mu_new[i],mu[i],reltol);
}
// saturated case
BOOST_CHECK_CLOSE(dmudp_diff,dmudp[0],reltol);
BOOST_CHECK_CLOSE(dmudr_diff,dmudr[0],reltol);
// unsaturated case
BOOST_CHECK_CLOSE(dmudp_diff_u,dmudp[3] , reltol);
BOOST_CHECK_CLOSE(dmudr_diff_u,dmudr[3] , reltol);
}
// test b
std::vector<double> b(n);
std::vector<double> B(n);
std::vector<double> invB(n);
std::vector<double> dinvBdp(n);
std::vector<double> dBdp(n);
std::vector<double> dbdr(n);
std::vector<double> dbdp(n);
double dbdp_diff;
double dbdr_diff;
double dbdp_diff_u;
double dbdr_diff_u;
for (int phase = 1; phase < 2; ++phase) {
props_[phase]->b(n, &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]);
//props_[phase]->B(n, p, z, B);
props_[phase]->dBdp(n, &p[0], &z[0], &B[0], &dBdp[0]);
dbdp_diff = (b[1]-b[0])/h_p;
dbdr_diff = (b[2]-b[0])/h_r;
dbdp_diff_u = (b[4]-b[3])/h_p;
dbdr_diff_u = (b[5]-b[3])/h_r;
for (int i = 0; i < n; ++i){
invB[i] = 1/B[i];
dinvBdp[i] = -1/pow(B[i],2) * dBdp[i];
}
for (int i = 0; i < n; ++i){
BOOST_CHECK_CLOSE(invB[i],b[i] , reltol);
BOOST_CHECK_CLOSE(dinvBdp[i],dbdp[i] , reltol);
}
// saturated case
BOOST_CHECK_CLOSE(dbdp_diff,dbdp[0], reltol);
BOOST_CHECK_CLOSE(dbdr_diff,dbdr[0], reltol);
// unsaturated case
BOOST_CHECK_CLOSE(dbdp_diff_u,dbdp[3], reltol);
BOOST_CHECK_CLOSE(dbdr_diff_u,dbdr[3], reltol);
}
// test bublepoint pressure
std::vector<double> rbub(n);
std::vector<double> drbubdp(n);
double drbubdp_diff;
double drbubdp_diff_u;
for (int phase = 1; phase < 2; ++phase) {
props_[phase] ->rbub(n, &p[0], &rbub[0], &drbubdp[0]);
drbubdp_diff = (rbub[1]-rbub[0])/h_p;
drbubdp_diff_u = (rbub[4]-rbub[3])/h_p;
// saturated case
BOOST_CHECK_CLOSE(drbubdp_diff,drbubdp[0], reltol);
// unsaturad case
BOOST_CHECK_CLOSE(drbubdp_diff_u,drbubdp[3], reltol);
}
testmu(reltol, n, np, p, r,z, props_, condition);
testb(reltol,n,np,p,r,z,props_,condition);
testrsSat(reltol,n,np,p,props_);
testrvSat(reltol,n,np,p,props_);
}
BOOST_AUTO_TEST_CASE(test_wetgas)
{
// read eclipse deck
const string filename = "wetgas.DATA";
cout << "Reading deck: " << filename << endl;
const EclipseGridParser deck (filename);
// setup pvt interface
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(deck,phase_usage_);
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
// approximation of the derivatives.
const int n = 6;
const int np = phase_usage_.num_phases;
// the tolerance for acceptable difference in values
const double reltol = 1e-9;
std::vector<double> p(n);
std::vector<double> r(n);
std::vector<PhasePresence> condition(n);
std::vector<double> z(n * np);
// Used for forward difference calculations
const double h_p = 1e4;
const double h_r = 1e-7;
// saturated
p[0] = 1e7;
p[1] = p[0] + h_p;
p[2] = p[0];
r[0] = 5e-5;
r[1] = 5e-5;
r[2] = 5e-5 + h_r;
condition[0].setFreeOil();
condition[1].setFreeOil();
condition[2].setFreeOil();
// undersaturated
p[3] = p[0];
p[4] = p[1];
p[5] = p[2];
r[3] = 1e-5;
r[4] = 1e-5;
r[5] = 1e-5 +h_r;
// Corresponing z factors, used to compare with the [p,z] interface
for (int i = 0; i < n; ++i) {
z[0+i*np] = 0; z[1+i*np] = r[i];
z[2+i*np] = 1;
}
testmu(reltol, n, np, p, r,z, props_, condition);
testb(reltol,n,np,p,r,z,props_,condition);
testrsSat(reltol,n,np,p,props_);
testrvSat(reltol,n,np,p,props_);
}
+71
View File
@@ -0,0 +1,71 @@
WATER
OIL
GAS
FIELD
-- PVT PROPERTIES OF DRY GAS (NO VAPOURISED OIL)
-- FROM SPE3 Blackoil Kleppe
--
-- 'Pressure' 'Oil FVF' 'Oil Visc'
PVDO
1214.7000 1.0632 0.3668
1814.7000 1.0518 0.4241
2414.7000 1.0418 0.5018
3014.7000 1.0332 0.6068
3214.7000 1.0308 0.6461
3364.7000 1.0291 0.6753
3414.7000 1.0285 0.6852
3443.8831 1.0282 0.6912
/
-- Wet Gas PVT Properties (Vapourised Oil)
-- Column Properties are:
-- 'Gas Pressure' 'Gas OGR' 'Gas FVF' 'Gas Visc'
-- Units: psia stb /Mscf rb /Mscf cp
PVTG
1214.7000 0.0013130 2.2799 0.0149
0 2.2815 0.01488/
1814.7000 0.00353 1.4401 0.01791
0.001313 1.4429 0.01782
0 1.4445 0.01735 /
2414.7000 0.01102 1.0438 0.02328
0.00353 1.0495 0.02267
0.001313 1.0512 0.0225
0 1.0522 0.02240 /
3014.7000 0.0331 0.8456 0.0318
0.01102 0.8489 0.02924
0.00353 0.8500 0.02844
0.001313 0.8503 0.02820
0 0.8505 0.02807 /
3214.7000 0.0454 0.8082 0.03539
0.0331 0.8080 0.03371
0.01102 0.8075 0.03113
0.00353 0.8073 0.03029
0.001313 0.8073 0.03004
0 0.8073 0.02989 /
3364.7000 0.05670 0.7875 0.0384
0.04540 0.7860 0.03667
0.03310 0.7843 0.03515
0.01102 0.7814 0.03429
0.00353 0.7804 0.03162
0.001313 0.7801 0.03136
0 0.7799 0.03121 /
3416.7575 0.0612 0.7816 0.03955
0.0567 0.7809 0.0386
0.0454 0.7789 0.03717
0.0331 0.7768 0.03564
0.01102 0.7731 0.03296
0.00353 0.7718 0.03207
0.001313 0.7714 0.03181
0 0.7712 0.03166 /
3449.3322 0.0642 0.7783 0.0403
0.0612 0.7777 0.0395
0.0567 0.7769 0.03892
0.0454 0.7747 0.03748
0.0331 0.7723 0.03594
0.01102 0.7681 0.03325
0.00353 0.7666 0.03236
0.001313 0.7662 0.0321
0 0.7660 0.03194 /
/