OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

allow temperature dependent PVT properties

Browse Source 
this has the nice side effect that non-uniform temperature fields also
become supported as long as they are "impressed" externally...
This commit is contained in:
Andreas Lauser
2014-11-20 18:57:29 +01:00
parent 7844931b41
commit d36d7e81cb
11 changed files with 74 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp
View File

@@ -127,6 +127,7 @@ namespace Opm {
struct SolutionState {
SolutionState(const int np);
ADB pressure;
ADB temperature;
std::vector<ADB> saturation;
ADB rs;
ADB rv;
@@ -277,6 +278,7 @@ namespace Opm {
ADB
fluidViscosity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@@ -285,6 +287,7 @@ namespace Opm {
ADB
fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@@ -293,6 +296,7 @@ namespace Opm {
ADB
fluidDensity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,

37
opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp
View File

@@ -365,6 +365,7 @@ namespace {
template<class T>
FullyImplicitBlackoilPolymerSolver<T>::SolutionState::SolutionState(const int np)
: pressure ( ADB::null())
, temperature( ADB::null())
, saturation(np, ADB::null())
, rs ( ADB::null())
, rv ( ADB::null())
@@ -422,6 +423,7 @@ namespace {
// automatically consistent with variableState() (and doing
// things automatically is all the rage in this module ;)
state.pressure = ADB::constant(state.pressure.value());
state.temperature = ADB::constant(state.temperature.value());
state.rs = ADB::constant(state.rs.value());
state.rv = ADB::constant(state.rv.value());
state.concentration = ADB::constant(state.concentration.value());
@@ -527,6 +529,10 @@ namespace {
int nextvar = 0;
state.pressure = vars[ nextvar++ ];
// Temperature.
const V temp = Eigen::Map<const V>(& x.temperature()[0], nc, 1);
state.pressure = ADB::constant(temp);
// Saturations
const std::vector<int>& bpat = vars[0].blockPattern();
{
@@ -592,6 +598,7 @@ namespace {
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB& press = state.pressure;
const ADB& temp = state.temperature;
const std::vector<ADB>& sat = state.saturation;
const ADB& rs = state.rs;
const ADB& rv = state.rv;
@@ -606,7 +613,7 @@ namespace {
for (int phase = 0; phase < maxnp; ++phase) {
if (active_[ phase ]) {
const int pos = pu.phase_pos[ phase ];
rq_[pos].b = fluidReciprocFVF(phase, pressure[phase], rs, rv, cond, cells_);
rq_[pos].b = fluidReciprocFVF(phase, pressure[phase], temp, rs, rv, cond, cells_);
rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos];
// DUMP(rq_[pos].b);
// DUMP(rq_[pos].accum[aix]);
@@ -666,6 +673,7 @@ namespace {
const std::vector<int> well_cells(wells_.well_cells, wells_.well_cells + nperf);
// Compute b, rsmax, rvmax values for perforations.
const ADB perf_press = subset(state.pressure, well_cells);
const ADB perf_temp = subset(state.temperature, well_cells);
std::vector<PhasePresence> perf_cond(nperf);
const std::vector<PhasePresence>& pc = phaseCondition();
for (int perf = 0; perf < nperf; ++perf) {
@@ -676,21 +684,21 @@ namespace {
std::vector<double> rssat_perf(nperf, 0.0);
std::vector<double> rvsat_perf(nperf, 0.0);
if (pu.phase_used[BlackoilPhases::Aqua]) {
const ADB bw = fluid_.bWat(perf_press, well_cells);
const ADB bw = fluid_.bWat(perf_press, perf_temp, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw.value();
}
assert(active_[Oil]);
const ADB perf_so = subset(state.saturation[pu.phase_pos[Oil]], well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) {
const ADB perf_rs = subset(state.rs, well_cells);
const ADB bo = fluid_.bOil(perf_press, perf_rs, perf_cond, well_cells);
const ADB bo = fluid_.bOil(perf_press, perf_temp, perf_rs, perf_cond, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo.value();
const V rssat = fluidRsSat(perf_press.value(), perf_so.value(), well_cells);
rssat_perf.assign(rssat.data(), rssat.data() + nperf);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = subset(state.rv, well_cells);
const ADB bg = fluid_.bGas(perf_press, perf_rv, perf_cond, well_cells);
const ADB bg = fluid_.bGas(perf_press, perf_temp, perf_rv, perf_cond, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg.value();
const V rvsat = fluidRvSat(perf_press.value(), perf_so.value(), well_cells);
rvsat_perf.assign(rvsat.data(), rvsat.data() + nperf);
@@ -1688,10 +1696,10 @@ namespace {
const std::vector<PhasePresence> cond = phaseCondition();
const ADB tr_mult = transMult(state.pressure);
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_);
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond, cells_);
rq_[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu;
const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_);
const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond, cells_);
ADB& head = rq_[phase].head;
@@ -1984,6 +1992,7 @@ namespace {
ADB
FullyImplicitBlackoilPolymerSolver<T>::fluidViscosity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@@ -1991,12 +2000,12 @@ namespace {
{
switch (phase) {
case Water:
return fluid_.muWat(p, cells);
return fluid_.muWat(p, temp, cells);
case Oil: {
return fluid_.muOil(p, rs, cond, cells);
return fluid_.muOil(p, temp, rs, cond, cells);
}
case Gas:
return fluid_.muGas(p, rv, cond, cells);
return fluid_.muGas(p, temp, rv, cond, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
@@ -2010,6 +2019,7 @@ namespace {
ADB
FullyImplicitBlackoilPolymerSolver<T>::fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@@ -2017,12 +2027,12 @@ namespace {
{
switch (phase) {
case Water:
return fluid_.bWat(p, cells);
return fluid_.bWat(p, temp, cells);
case Oil: {
return fluid_.bOil(p, rs, cond, cells);
return fluid_.bOil(p, temp, rs, cond, cells);
}
case Gas:
return fluid_.bGas(p, rv, cond, cells);
return fluid_.bGas(p, temp, rv, cond, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
@@ -2036,13 +2046,14 @@ namespace {
ADB
FullyImplicitBlackoilPolymerSolver<T>::fluidDensity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
const double* rhos = fluid_.surfaceDensity();
ADB b = fluidReciprocFVF(phase, p, rs, rv, cond, cells);
ADB b = fluidReciprocFVF(phase, p, temp, rs, rv, cond, cells);
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
if (phase == Oil && active_[Gas]) {
// It is correct to index into rhos with canonical phase indices.

33
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
View File

@@ -263,6 +263,7 @@ namespace {
FullyImplicitCompressiblePolymerSolver::SolutionState::SolutionState(const int np)
: pressure ( ADB::null())
, temperature( ADB::null())
, saturation(np, ADB::null())
, concentration( ADB::null())
, qs ( ADB::null())
@@ -331,6 +332,10 @@ namespace {
const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
state.pressure = ADB::constant(p, bpat);
// Temperature.
const V T = Eigen::Map<const V>(& x.temperature()[0], nc, 1);
state.temperature = ADB::constant(T);
// Saturation.
assert (not x.saturation().empty());
const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
@@ -449,6 +454,7 @@ namespace {
{
const ADB& press = state.pressure;
const ADB& temp = state.temperature;
const std::vector<ADB>& sat = state.saturation;
const ADB& c = state.concentration;
@@ -458,7 +464,7 @@ namespace {
const ADB pv_mult = poroMult(press);
for (int phase = 0; phase < 2; ++phase) {
rq_[phase].b = fluidReciprocFVF(phase, pressure[phase], cond, cells_);
rq_[phase].b = fluidReciprocFVF(phase, pressure[phase], temp, cond, cells_);
}
rq_[0].accum[aix] = pv_mult * rq_[0].b * sat[0];
rq_[1].accum[aix] = pv_mult * rq_[1].b * sat[1];
@@ -566,16 +572,17 @@ namespace {
}
ADB cell_rho_total = ADB::constant(V::Zero(nc), state.pressure.blockPattern());
std::vector<ADB> press = computePressures(state);
const ADB& temp = state.temperature;
const std::vector<PhasePresence> cond = phaseCondition();
for (int phase = 0; phase < 2; ++phase) {
const ADB cell_rho = fluidDensity(phase, press[phase], cond, cells_);
const ADB cell_rho = fluidDensity(phase, press[phase], temp, cond, cells_);
cell_rho_total += state.saturation[phase] * cell_rho;
}
ADB inj_rho_total = ADB::constant(V::Zero(nperf), state.pressure.blockPattern());
assert(np == wells_.number_of_phases);
const DataBlock compi = Eigen::Map<const DataBlock>(wells_.comp_frac, nw, np);
for (int phase = 0; phase < 2; ++phase) {
const ADB cell_rho = fluidDensity(phase, press[phase], cond, cells_);
const ADB cell_rho = fluidDensity(phase, press[phase], temp, cond, cells_);
const V fraction = compi.col(phase);
inj_rho_total += (wops_.w2p * fraction.matrix()).array() * subset(cell_rho, well_cells);
}
@@ -840,15 +847,16 @@ namespace {
const ADB tr_mult = transMult(state.pressure);
const std::vector<PhasePresence> cond = phaseCondition();
std::vector<ADB> press = computePressures(state);
const ADB& temp = state.temperature;
const ADB mu_w = fluidViscosity(0, press[0], cond, cells_);
const ADB mu_w = fluidViscosity(0, press[0], temp, cond, cells_);
ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu_w.value().data());
rq_[0].mob = tr_mult * krw_eff * inv_wat_eff_vis;
rq_[2].mob = tr_mult * mc * krw_eff * inv_wat_eff_vis;
const ADB mu_o = fluidViscosity(1, press[1], cond, cells_);
const ADB mu_o = fluidViscosity(1, press[1], temp, cond, cells_);
rq_[1].mob = tr_mult * kro / mu_o;
for (int phase = 0; phase < 2; ++phase) {
const ADB rho = fluidDensity(phase, press[phase], cond, cells_);
const ADB rho = fluidDensity(phase, press[phase], temp, cond, cells_);
ADB& head = rq_[ phase ].head;
// compute gravity potensial using the face average as in eclipse and MRST
const ADB rhoavg = ops_.caver * rho;
@@ -893,15 +901,16 @@ namespace {
ADB
FullyImplicitCompressiblePolymerSolver::fluidViscosity(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
const ADB null = ADB::constant(V::Zero(grid_.number_of_cells, 1), p.blockPattern());
switch (phase) {
case Water:
return fluid_.muWat(p, cells);
return fluid_.muWat(p, T, cells);
case Oil: {
return fluid_.muOil(p, null, cond, cells);
return fluid_.muOil(p, T, null, cond, cells);
}
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
@@ -915,15 +924,16 @@ namespace {
ADB
FullyImplicitCompressiblePolymerSolver::fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
const ADB null = ADB::constant(V::Zero(grid_.number_of_cells, 1), p.blockPattern());
switch (phase) {
case Water:
return fluid_.bWat(p, cells);
return fluid_.bWat(p, T, cells);
case Oil: {
return fluid_.bOil(p, null, cond, cells);
return fluid_.bOil(p, T, null, cond, cells);
}
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
@@ -937,11 +947,12 @@ namespace {
ADB
FullyImplicitCompressiblePolymerSolver::fluidDensity(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
const double* rhos = fluid_.surfaceDensity();
ADB b = fluidReciprocFVF(phase, p, cond, cells);
ADB b = fluidReciprocFVF(phase, p, T, cond, cells);
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
return rho;
}

4
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp
View File

@@ -107,6 +107,7 @@ namespace Opm {
struct SolutionState {
SolutionState(const int np);
ADB pressure;
ADB temperature;
std::vector<ADB> saturation;
ADB concentration;
ADB qs;
@@ -213,18 +214,21 @@ namespace Opm {
ADB
fluidViscosity(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const;
ADB
fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const;
ADB
fluidDensity(const int phase,
const ADB& p ,
const ADB& T ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const;