diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 8d6ce4376..9eceb3565 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -454,6 +454,8 @@ namespace Opm { computeMassFlux(const int actph , const V& transi, const ADB& kr , + const ADB& mu , + const ADB& rho , const ADB& p , const SolutionState& state ); diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 0164461f9..d290b807a 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -752,7 +752,7 @@ namespace detail { for (int phase = 0; phase < maxnp; ++phase) { if (active_[ phase ]) { const int pos = pu.phase_pos[ phase ]; - rq_[pos].b = fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond); + rq_[pos].b = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond); rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos]; // OPM_AD_DUMP(rq_[pos].b); // OPM_AD_DUMP(rq_[pos].accum[aix]); @@ -989,7 +989,10 @@ namespace detail { const std::vector kr = asImpl().computeRelPerm(state); #pragma omp parallel for schedule(static) for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) { - asImpl().computeMassFlux(phaseIdx, trans_all, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state); + const std::vector& cond = phaseCondition(); + const ADB mu = asImpl().fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond); + const ADB rho = asImpl().fluidDensity(canph_[phaseIdx], rq_[phaseIdx].b, state.rs, state.rv); + asImpl().computeMassFlux(phaseIdx, trans_all, kr[canph_[phaseIdx]], mu, rho, state.canonical_phase_pressures[canph_[phaseIdx]], state); residual_.material_balance_eq[ phaseIdx ] = pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0]) @@ -2404,25 +2407,21 @@ namespace detail { - - template void BlackoilModelBase::computeMassFlux(const int actph , const V& transi, const ADB& kr , + const ADB& mu , + const ADB& rho , const ADB& phasePressure, const SolutionState& state) { // Compute and store mobilities. - const int canonicalPhaseIdx = canph_[ actph ]; - const std::vector& cond = phaseCondition(); const ADB tr_mult = transMult(state.pressure); - const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv, cond); rq_[ actph ].mob = tr_mult * kr / mu; // Compute head differentials. Gravity potential is done using the face average as in eclipse and MRST. - const ADB rho = fluidDensity(canonicalPhaseIdx, rq_[actph].b, state.rs, state.rv); const ADB rhoavg = ops_.caver * rho; rq_[ actph ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix())); if (use_threshold_pressure_) { diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index 4e0fe7166..e8cc9a063 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -63,6 +63,7 @@ namespace Opm { /// \param[in] has_vapoil turn on vaporized oil feature /// \param[in] terminal_output request output to cout/cerr /// \param[in] has_solvent turn on solvent feature + /// \param[in] is_miscible turn on miscible feature BlackoilSolventModel(const typename Base::ModelParameters& param, const Grid& grid, const BlackoilPropsAdInterface& fluid, @@ -75,7 +76,8 @@ namespace Opm { const bool has_disgas, const bool has_vapoil, const bool terminal_output, - const bool has_solvent); + const bool has_solvent, + const bool is_miscible); /// Apply an update to the primary variables, chopped if appropriate. /// \param[in] dx updates to apply to primary variables @@ -106,6 +108,10 @@ namespace Opm { const bool has_solvent_; const int solvent_pos_; const SolventPropsAdFromDeck& solvent_props_; + const bool is_miscible_; + std::vector mu_eff_; + std::vector b_eff_; + // Need to declare Base members we want to use here. using Base::grid_; @@ -141,9 +147,6 @@ namespace Opm { using Base::computePressures; using Base::computeGasPressure; using Base::applyThresholdPressures; - using Base::fluidViscosity; - using Base::fluidReciprocFVF; - using Base::fluidDensity; using Base::fluidRsSat; using Base::fluidRvSat; using Base::poroMult; @@ -162,6 +165,28 @@ namespace Opm { std::vector computeRelPerm(const SolutionState& state) const; + ADB + fluidViscosity(const int phase, + const ADB& p , + const ADB& temp , + const ADB& rs , + const ADB& rv , + const std::vector& cond) const; + + ADB + fluidReciprocFVF(const int phase, + const ADB& p , + const ADB& temp , + const ADB& rs , + const ADB& rv , + const std::vector& cond) const; + + ADB + fluidDensity(const int phase, + const ADB& b, + const ADB& rs, + const ADB& rv) const; + void makeConstantState(SolutionState& state) const; @@ -198,12 +223,19 @@ namespace Opm { computeMassFlux(const int actph , const V& transi, const ADB& kr , + const ADB& mu , + const ADB& rho , const ADB& p , const SolutionState& state ); const std::vector phaseCondition() const {return this->phaseCondition_;} + // compute effective viscosities (mu_eff_) and effective b factors (b_eff_) using the ToddLongstaff model + void computeEffectiveProperties(const SolutionState& state); + + // compute density and viscosity using the ToddLongstaff mixing model + void computeToddLongstaffMixing(std::vector& viscosity, std::vector& density, const std::vector& saturations, const Opm::PhaseUsage pu); }; diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index 79e848d35..9a2c5fee5 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -82,12 +82,15 @@ namespace Opm { const bool has_disgas, const bool has_vapoil, const bool terminal_output, - const bool has_solvent) + const bool has_solvent, + const bool is_miscible) : Base(param, grid, fluid, geo, rock_comp_props, wells_arg, linsolver, eclState, has_disgas, has_vapoil, terminal_output), has_solvent_(has_solvent), solvent_pos_(detail::solventPos(fluid.phaseUsage())), - solvent_props_(solvent_props) + solvent_props_(solvent_props), + is_miscible_(is_miscible) + { if (has_solvent_) { @@ -101,7 +104,10 @@ namespace Opm { } residual_.matbalscale.resize(fluid_.numPhases() + 1, 0.0031); // use the same as gas - + } + if (is_miscible_) { + mu_eff_.resize(fluid_.numPhases() + 1, ADB::null()); + b_eff_.resize(fluid_.numPhases() + 1, ADB::null()); } } @@ -203,7 +209,8 @@ namespace Opm { const Opm::PhaseUsage& pu = fluid_.phaseUsage(); const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[Gas]]; - rq_[solvent_pos_].b = solvent_props_.bSolvent(pg,cells_); + const std::vector& cond = phaseCondition(); + rq_[solvent_pos_].b = fluidReciprocFVF(Solvent, pg, state.temperature, state.rs, state.rv,cond); rq_[solvent_pos_].accum[aix] = pv_mult * rq_[solvent_pos_].b * ss; } } @@ -354,33 +361,42 @@ namespace Opm { DataBlock b(nperf, pu.num_phases); std::vector rsmax_perf(nperf, 0.0); std::vector rvmax_perf(nperf, 0.0); + const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); if (pu.phase_used[BlackoilPhases::Aqua]) { - const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; } + assert(active_[Oil]); + assert(active_[Gas]); + const ADB perf_rv = subset(state.rv, well_cells); + const ADB perf_rs = subset(state.rs, well_cells); const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); if (pu.phase_used[BlackoilPhases::Liquid]) { - const ADB perf_rs = subset(state.rs, well_cells); const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); + //const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value(); b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; const V rssat = fluidRsSat(avg_press, perf_so, well_cells); rsmax_perf.assign(rssat.data(), rssat.data() + nperf); } V surf_dens_copy = superset(fluid_.surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); for (int phase = 1; phase < pu.num_phases; ++phase) { - if ( phase != pu.phase_pos[BlackoilPhases::Vapour]) { + if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) { continue; // the gas surface density is added after the solvent is accounted for. } surf_dens_copy += superset(fluid_.surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); } if (pu.phase_used[BlackoilPhases::Vapour]) { - const ADB perf_rv = subset(state.rv, well_cells); + // Unclear wether the effective or the pure values should be used for the wells + // the current usage of unmodified properties values gives best match. + //V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value(); V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); V rhog = fluid_.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells); if (has_solvent_) { + const V bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value(); + //const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value(); + // A weighted sum of the b-factors of gas and solvent are used. const int nc = Opm::AutoDiffGrid::numCells(grid_); @@ -404,6 +420,7 @@ namespace Opm { } } } + F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction; bg = bg * (ones - F_solvent); @@ -510,53 +527,114 @@ namespace Opm { BlackoilSolventModel::computeMassFlux(const int actph , const V& transi, const ADB& kr , + const ADB& mu , + const ADB& rho , const ADB& phasePressure, const SolutionState& state) { - Base::computeMassFlux(actph, transi, kr, phasePressure, state); const int canonicalPhaseIdx = canph_[ actph ]; + // make a copy to make it possible to modify it + ADB kr_mod = kr; if (canonicalPhaseIdx == Gas) { if (has_solvent_) { const int nc = Opm::UgGridHelpers::numCells(grid_); - const Opm::PhaseUsage& pu = fluid_.phaseUsage(); const ADB zero = ADB::constant(V::Zero(nc)); + const V ones = V::Constant(nc, 1.0); + const ADB& ss = state.solvent_saturation; const ADB& sg = (active_[ Gas ] ? state.saturation[ pu.phase_pos[ Gas ] ] : zero); - Selector zero_selector(ss.value() + sg.value(), Selector::Zero); - ADB F_solvent = zero_selector.select(ss, ss / (ss + sg)); - V ones = V::Constant(nc, 1.0); + const ADB F_solvent = zero_selector.select(zero, ss / (ss + sg)); - const ADB tr_mult = transMult(state.pressure); - const ADB mu = solvent_props_.muSolvent(phasePressure,cells_); + // Compute solvent properties + const std::vector& cond = phaseCondition(); + ADB mu_s = fluidViscosity(Solvent, phasePressure,state.temperature, state.rs, state.rv, cond); + ADB rho_s = fluidDensity(Solvent,rq_[solvent_pos_].b, state.rs, state.rv); - rq_[solvent_pos_].mob = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * tr_mult * kr / mu; + // Compute solvent relperm and mass flux + ADB krs = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * kr_mod; + Base::computeMassFlux(solvent_pos_, transi, krs, mu_s, rho_s, phasePressure, state); - const ADB rho_solvent = solvent_props_.solventSurfaceDensity(cells_) * rq_[solvent_pos_].b; - const ADB rhoavg_solvent = ops_.caver * rho_solvent; - rq_[ solvent_pos_ ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg_solvent * (ops_.ngrad * geo_.z().matrix())); + // Modify gas relperm + kr_mod = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * kr_mod; - UpwindSelector upwind_solvent(grid_, ops_, rq_[solvent_pos_].dh.value()); - // Compute solvent flux. - rq_[solvent_pos_].mflux = upwind_solvent.select(rq_[solvent_pos_].b * rq_[solvent_pos_].mob) * (transi * rq_[solvent_pos_].dh); - - // Update gas mobility and flux - rq_[actph].mob = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * rq_[actph].mob; - - const ADB& b = rq_[ actph ].b; - const ADB& mob = rq_[ actph ].mob; - const ADB& dh = rq_[ actph ].dh; - UpwindSelector upwind_gas(grid_, ops_, dh.value()); - rq_[ actph ].mflux = upwind_gas.select(b * mob) * (transi * dh); } } + // Compute mobility and flux + Base::computeMassFlux(actph, transi, kr_mod, mu, rho, phasePressure, state); } + template + ADB + BlackoilSolventModel::fluidViscosity(const int phase, + const ADB& p , + const ADB& temp , + const ADB& rs , + const ADB& rv , + const std::vector& cond) const + { + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + if (phase == Solvent) { + if (!is_miscible_) { + return solvent_props_.muSolvent(p, cells_); + } else { + return mu_eff_[solvent_pos_]; + } + + } else { + if (!is_miscible_) { + return Base::fluidViscosity(phase, p, temp, rs, rv, cond); + } else { + return mu_eff_[pu.phase_pos[ phase ]]; + } + } + } + + template + ADB + BlackoilSolventModel::fluidReciprocFVF(const int phase, + const ADB& p , + const ADB& temp , + const ADB& rs , + const ADB& rv , + const std::vector& cond) const + { + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + if (phase == Solvent) { + if (!is_miscible_) { + return solvent_props_.bSolvent(p, cells_); + } else { + return b_eff_[solvent_pos_]; + } + + } else { + if (!is_miscible_) { + return Base::fluidReciprocFVF(phase, p, temp, rs, rv, cond); + } else { + return b_eff_[pu.phase_pos[ phase ]]; + } + } + } + + template + ADB + BlackoilSolventModel::fluidDensity(const int phase, + const ADB& b, + const ADB& rs, + const ADB& rv) const + { + if (phase == Solvent && has_solvent_) { + return solvent_props_.solventSurfaceDensity(cells_) * b; + } else { + return Base::fluidDensity(phase, b, rs, rv); + } + } + template std::vector BlackoilSolventModel::computeRelPerm(const SolutionState& state) const @@ -581,13 +659,227 @@ namespace Opm { if (has_solvent_) { const ADB& ss = state.solvent_saturation; - return fluid_.relperm(sw, so, sg+ss, cells_); + if (is_miscible_) { + + std::vector relperm = fluid_.relperm(sw, so, sg+ss, cells_); + + Selector zero_selector(ss.value() + sg.value(), Selector::Zero); + ADB F_solvent = zero_selector.select(ss, ss / (ss + sg)); + const ADB misc = solvent_props_.miscibilityFunction(F_solvent, cells_); + + assert(active_[ Oil ]); + assert(active_[ Gas ]); + + const ADB sn = ss + so + sg; + + // adjust endpoints + const V sgcr = fluid_.scaledCriticalGasSaturations(cells_); + const V sogcr = fluid_.scaledCriticalOilinGasSaturations(cells_); + const ADB sorwmis = solvent_props_.miscibleResidualOilSaturationFunction(sw, cells_); + const ADB sgcwmis = solvent_props_.miscibleCriticalGasSaturationFunction(sw, cells_); + + const V ones = V::Constant(nc, 1.0); + ADB sor = misc * sorwmis + (ones - misc) * sogcr; + ADB sgc = misc * sgcwmis + (ones - misc) * sgcr; + + // avoid negative values + Selector negative_selector(sgc.value(), Selector::LessEqualZero); + sgc = negative_selector.select(zero, sgc); + + const ADB ssg = ss + sg - sgc; + const ADB sn_eff = sn - sor - sgc; + + // avoid negative values + Selector zeroSn_selector(sn_eff.value(), Selector::Zero); + const ADB F_totalGas = zeroSn_selector.select( zero, ssg / sn_eff); + + const ADB mkrgt = solvent_props_.miscibleSolventGasRelPermMultiplier(F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_); + const ADB mkro = solvent_props_.miscibleOilRelPermMultiplier(ones - F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_); + + const V eps = V::Constant(nc, 1e-7); + Selector noOil_selector(so.value()-eps, Selector::LessEqualZero); + relperm[Gas] = noOil_selector.select(relperm[Gas], (ones - misc) * relperm[Gas] + misc * mkrgt); + relperm[Oil] = noOil_selector.select(relperm[Oil], (ones - misc) * relperm[Oil] + misc * mkro); + return relperm; + } else { + return fluid_.relperm(sw, so, sg+ss, cells_); + } } else { return fluid_.relperm(sw, so, sg, cells_); } } + template + void + BlackoilSolventModel::computeEffectiveProperties(const SolutionState& state) + { + // Viscosity + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + const int np = fluid_.numPhases(); + const int nc = Opm::UgGridHelpers::numCells(grid_); + const ADB zero = ADB::constant(V::Zero(nc)); + + const ADB& pw = state.canonical_phase_pressures[pu.phase_pos[Water]]; + const ADB& po = state.canonical_phase_pressures[pu.phase_pos[Oil]]; + const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[Gas]]; + const std::vector& cond = phaseCondition(); + + const ADB mu_w = fluid_.muWat(pw, state.temperature, cells_); + const ADB mu_o = fluid_.muOil(po, state.temperature, state.rs, cond, cells_); + const ADB mu_g = fluid_.muGas(pg, state.temperature, state.rv, cond, cells_); + const ADB mu_s = solvent_props_.muSolvent(pg,cells_); + std::vector viscosity(np + 1, ADB::null()); + viscosity[pu.phase_pos[Oil]] = mu_o; + viscosity[pu.phase_pos[Gas]] = mu_g; + viscosity[pu.phase_pos[Water]] = mu_w; + viscosity[solvent_pos_] = mu_s; + + // Density + const ADB bw = fluid_.bWat(pw, state.temperature, cells_); + const ADB bo = fluid_.bOil(po, state.temperature, state.rs, cond, cells_); + const ADB bg = fluid_.bGas(pg, state.temperature, state.rv, cond, cells_); + const ADB bs = solvent_props_.bSolvent(pg, cells_); + + std::vector density(np + 1, ADB::null()); + density[pu.phase_pos[Oil]] = fluidDensity(Oil, bo, state.rs, state.rv); + density[pu.phase_pos[Gas]] = fluidDensity(Gas, bg, state.rs, state.rv); + density[pu.phase_pos[Water]] = fluidDensity(Water, bw, state.rs, state.rv); + density[solvent_pos_] = fluidDensity(Solvent, bs, state.rs, state.rv); + + // Effective saturations + const ADB& ss = state.solvent_saturation; + const ADB& so = state.saturation[ pu.phase_pos[ Oil ] ]; + const ADB& sg = (active_[ Gas ] + ? state.saturation[ pu.phase_pos[ Gas ] ] + : zero); + const ADB& sw = (active_[ Water ] + ? state.saturation[ pu.phase_pos[ Water ] ] + : zero); + + const ADB sorwmis = solvent_props_.miscibleResidualOilSaturationFunction(sw, cells_); + const ADB sgcwmis = solvent_props_.miscibleCriticalGasSaturationFunction(sw, cells_); + + std::vector effective_saturations (np + 1, ADB::null()); + effective_saturations[pu.phase_pos[Oil]] = so - sorwmis; + effective_saturations[pu.phase_pos[Gas]] = sg - sgcwmis; + effective_saturations[pu.phase_pos[Water]] = sw; + effective_saturations[solvent_pos_] = ss - sgcwmis; + + // Compute effective viscosities and densities + computeToddLongstaffMixing(viscosity, density, effective_saturations, pu); + + // Store the computed volume factors and viscosities + b_eff_[pu.phase_pos[ Water ]] = bw; + b_eff_[pu.phase_pos[ Oil ]] = density[pu.phase_pos[ Oil ]] / (fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * state.rs); + b_eff_[pu.phase_pos[ Gas ]] = density[pu.phase_pos[ Gas ]] / (fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * state.rv); + b_eff_[solvent_pos_] = density[solvent_pos_] / solvent_props_.solventSurfaceDensity(cells_); + + mu_eff_[pu.phase_pos[ Water ]] = mu_w; + mu_eff_[pu.phase_pos[ Oil ]] = viscosity[pu.phase_pos[ Oil ]]; + mu_eff_[pu.phase_pos[ Gas ]] = viscosity[pu.phase_pos[ Gas ]]; + mu_eff_[solvent_pos_] = viscosity[solvent_pos_]; + } + + template + void + BlackoilSolventModel::computeToddLongstaffMixing(std::vector& viscosity, std::vector& density, const std::vector& saturations, const Opm::PhaseUsage pu) + { + const int nc = cells_.size(); + const V ones = V::Constant(nc, 1.0); + const ADB zero = ADB::constant(V::Zero(nc)); + + // Calculation of effective saturations + ADB so_eff = saturations[pu.phase_pos[ Oil ]]; + ADB sg_eff = saturations[pu.phase_pos[ Gas ]]; + ADB ss_eff = saturations[solvent_pos_]; + + // Avoid negative values + Selector negative_selectorSo(so_eff.value(), Selector::LessZero); + Selector negative_selectorSg(sg_eff.value(), Selector::LessZero); + Selector negative_selectorSs(ss_eff.value(), Selector::LessZero); + so_eff = negative_selectorSo.select(zero, so_eff); + sg_eff = negative_selectorSg.select(zero, sg_eff); + ss_eff = negative_selectorSs.select(zero, ss_eff); + + // Make copy of the pure viscosities + const ADB mu_o = viscosity[pu.phase_pos[ Oil ]]; + const ADB mu_g = viscosity[pu.phase_pos[ Gas ]]; + const ADB mu_s = viscosity[solvent_pos_]; + + const ADB sn_eff = so_eff + sg_eff + ss_eff; + const ADB sos_eff = so_eff + ss_eff; + const ADB ssg_eff = ss_eff + sg_eff; + + // Avoid division by zero + Selector zero_selectorSos(sos_eff.value(), Selector::Zero); + Selector zero_selectorSsg(ssg_eff.value(), Selector::Zero); + Selector zero_selectorSn(sn_eff.value(), Selector::Zero); + + const ADB mu_s_pow = pow(mu_s,0.25); + const ADB mu_o_pow = pow(mu_o,0.25); + const ADB mu_g_pow = pow(mu_g,0.25); + + const ADB mu_mos = zero_selectorSos.select(mu_o + mu_s, mu_o * mu_s / pow( ( (so_eff / sos_eff) * mu_s_pow) + ( (ss_eff / sos_eff) * mu_o_pow) , 4.0)); + const ADB mu_msg = zero_selectorSsg.select(mu_g + mu_s , mu_g * mu_s / pow( ( (sg_eff / ssg_eff) * mu_s_pow) + ( (ss_eff / ssg_eff) * mu_g_pow) , 4.0)); + const ADB mu_m = zero_selectorSn.select(mu_s + mu_o + mu_g, mu_o * mu_s * mu_g / pow( ( (so_eff / sn_eff) * mu_s_pow * mu_g_pow) + + ( (ss_eff / sn_eff) * mu_o_pow * mu_g_pow) + ( (sg_eff / sn_eff) * mu_s_pow * mu_o_pow), 4.0)); + // Mixing parameter for viscosity + const double mix_param_mu = solvent_props_.mixingParameterViscosity(); + + // Update viscosities + viscosity[pu.phase_pos[ Oil ]] = pow(mu_o,1.0 - mix_param_mu) * pow(mu_mos,mix_param_mu); + viscosity[pu.phase_pos[ Gas ]] = pow(mu_g,1.0 - mix_param_mu) * pow(mu_msg,mix_param_mu); + viscosity[solvent_pos_] = pow(mu_s,1.0 - mix_param_mu) * pow(mu_m,mix_param_mu); + + // Density + ADB& rho_o = density[pu.phase_pos[ Oil ]]; + ADB& rho_g = density[pu.phase_pos[ Gas ]]; + ADB& rho_s = density[solvent_pos_]; + + // mixing parameter for density + const double mix_param_rho = solvent_props_.mixingParameterDensity(); + + // compute effective viscosities for density calculations. These have to + // be recomputed as a different mixing parameter may be used. + const ADB mu_o_eff = pow(mu_o,1.0 - mix_param_rho) * pow(mu_mos,mix_param_rho); + const ADB mu_g_eff = pow(mu_g,1.0 - mix_param_rho) * pow(mu_msg,mix_param_rho); + const ADB mu_s_eff = pow(mu_s,1.0 - mix_param_rho) * pow(mu_m,mix_param_rho); + + const ADB sog_eff = so_eff + sg_eff; + const ADB sof = so_eff / sog_eff; + const ADB sgf = sg_eff / sog_eff; + + // Effective densities + const ADB mu_sog_pow = mu_s_pow * ( (sgf * mu_o_pow) + (sof * mu_g_pow) ); + const ADB mu_o_eff_pow = pow(mu_o_eff,0.25); + const ADB mu_g_eff_pow = pow(mu_g_eff,0.25); + const ADB mu_s_eff_pow = pow(mu_s_eff,0.25); + const ADB sfraction_oe = (mu_o_pow * (mu_o_eff_pow - mu_s_pow)) / (mu_o_eff_pow * (mu_o_pow - mu_s_pow)); + const ADB sfraction_ge = (mu_g_pow * (mu_s_pow - mu_g_eff_pow)) / (mu_g_eff_pow * (mu_s_pow - mu_g_pow)); + const ADB sfraction_se = (mu_sog_pow - ( mu_o_pow * mu_g_pow * mu_s_pow / mu_s_eff_pow) ) / ( mu_sog_pow - (mu_o_pow * mu_g_pow)); + const ADB rho_o_eff = (rho_o * sfraction_oe) + (rho_s * (ones - sfraction_oe)); + const ADB rho_g_eff = (rho_g * sfraction_ge) + (rho_s * (ones - sfraction_ge)); + const ADB rho_s_eff = (rho_s * sfraction_se) + (rho_g * sgf * (ones - sfraction_se)) + (rho_o * sof * (ones - sfraction_se)); + + // Avoid division by zero for equal mobilities. For equal mobilities the effecitive density is calculated + // based on the saturation fraction directly. + Selector unitGasSolventMobilityRatio_selector(mu_s.value() - mu_g.value(), Selector::Zero); + Selector unitOilSolventMobilityRatio_selector(mu_s.value() - mu_o.value(), Selector::Zero); + + // Effective densities when the mobilities are equal + const ADB rho_m = zero_selectorSn.select(zero, (rho_o * so_eff / sn_eff) + (rho_g * sg_eff / sn_eff) + (rho_s * ss_eff / sn_eff)); + const ADB rho_o_eff_simple = ((ones - mix_param_rho) * rho_o) + (mix_param_rho * rho_m); + const ADB rho_g_eff_simple = ((ones - mix_param_rho) * rho_g) + (mix_param_rho * rho_m); + //const ADB rho_s_eff_simple = ((ones - mix_param_rho) * rho_s) + (mix_param_rho * rho_m); + + // Update densities + rho_o = unitOilSolventMobilityRatio_selector.select(rho_o_eff_simple, rho_o_eff); + rho_g = unitGasSolventMobilityRatio_selector.select(rho_g_eff_simple, rho_g_eff); + rho_s = rho_s_eff; + + } + template void @@ -611,9 +903,16 @@ namespace Opm { makeConstantState(state0); // Compute initial accumulation contributions // and well connection pressures. + if (is_miscible_) { + computeEffectiveProperties(state0); + } + computeAccum(state0, 0); computeWellConnectionPressures(state0, well_state); } + if (is_miscible_) { + computeEffectiveProperties(state); + } // -------- Mass balance equations -------- assembleMassBalanceEq(state); @@ -661,14 +960,15 @@ namespace Opm { Selector zero_selector(ss.value() + sg.value(), Selector::Zero); ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells); V ones = V::Constant(nperf,1.0); + b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos]; b_perfcells[gas_pos] += (F_solvent * subset(rq_[solvent_pos_].b, well_cells)); + } if (param_.solve_welleq_initially_ && initial_assembly) { // solve the well equations as a pre-processing step Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state); } - Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state); Base::addWellFluxEq(cq_s, state); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp index 33802296b..bad80a289 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp @@ -131,6 +131,7 @@ namespace Opm bool has_solvent_; DeckConstPtr deck_; SolventPropsAdFromDeck solvent_props_; + bool is_miscible_; }; diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp index 3a3622acc..694471752 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp @@ -54,9 +54,10 @@ namespace Opm , has_solvent_(has_solvent) , deck_(deck) , solvent_props_(solvent_props) + , is_miscible_(false) { if(deck->hasKeyword("MISCIBLE")) { - std::cerr << "MISICIBLE keyword is present. Mixing is not currently supported" << std::endl; + is_miscible_ = true; } } @@ -80,7 +81,8 @@ namespace Opm BaseType::has_disgas_, BaseType::has_vapoil_, BaseType::terminal_output_, - has_solvent_)); + has_solvent_, + is_miscible_)); if (!BaseType::threshold_pressures_by_face_.empty()) { model->setThresholdPressures(BaseType::threshold_pressures_by_face_); diff --git a/opm/autodiff/SolventPropsAdFromDeck.cpp b/opm/autodiff/SolventPropsAdFromDeck.cpp index ef107aa71..7bcb2f824 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.cpp +++ b/opm/autodiff/SolventPropsAdFromDeck.cpp @@ -23,6 +23,7 @@ #include #include +#include namespace Opm @@ -96,7 +97,7 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, // relative permeabilty multiplier if (!ssfnTables.empty()) { - int numRegions = pvdsTables.size(); + int numRegions = ssfnTables.size(); if(numRegions > 1) { OPM_THROW(std::runtime_error, "Only single table saturation function supported for SSFN"); @@ -119,7 +120,157 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, } else { OPM_THROW(std::runtime_error, "SSFN must be specified in SOLVENT runs\n"); } + + + if (deck->hasKeyword("MISCIBLE") ) { + // misicible hydrocabon relative permeability wrt water + const TableContainer& sof2Tables = tables->getSof2Tables(); + if (!sof2Tables.empty()) { + + int numRegions = sof2Tables.size(); + + if(numRegions > 1) { + OPM_THROW(std::runtime_error, "Only single table saturation function supported for SOF2"); + } + // resize the attributes of the object + krn_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + const Opm::Sof2Table& sof2Table = sof2Tables.getTable(regionIdx); + + // Copy data + // Sn = So + Sg + Ss; + const auto& sn = sof2Table.getSoColumn(); + const auto& krn = sof2Table.getKroColumn(); + + krn_[regionIdx] = NonuniformTableLinear(sn, krn); + } + + } else { + OPM_THROW(std::runtime_error, "SOF2 must be specified in MISCIBLE (SOLVENT) runs\n"); + } + + const TableContainer& miscTables = tables->getMiscTables(); + if (!miscTables.empty()) { + + int numRegions = miscTables.size(); + + if(numRegions > 1) { + OPM_THROW(std::runtime_error, "Only single table miscibility function supported for MISC"); + } + // resize the attributes of the object + misc_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + const Opm::MiscTable& miscTable = miscTables.getTable(regionIdx); + + // Copy data + // solventFraction = Ss / (Ss + Sg); + const auto& solventFraction = miscTable.getSolventFractionColumn(); + const auto& misc = miscTable.getMiscibilityColumn(); + + misc_[regionIdx] = NonuniformTableLinear(solventFraction, misc); + + } + } else { + OPM_THROW(std::runtime_error, "MISC must be specified in MISCIBLE (SOLVENT) runs\n"); + } + + // miscible relative permeability multipleiers + const TableContainer& msfnTables = tables->getMsfnTables(); + if (!msfnTables.empty()) { + + int numRegions = msfnTables.size(); + + if(numRegions > 1) { + OPM_THROW(std::runtime_error, "Only single table saturation function supported for MSFN"); + } + // resize the attributes of the object + mkrsg_.resize(numRegions); + mkro_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + const Opm::MsfnTable& msfnTable = msfnTables.getTable(regionIdx); + + // Copy data + // Ssg = Ss + Sg; + const auto& Ssg = msfnTable.getGasPhaseFractionColumn(); + const auto& krsg = msfnTable.getGasSolventRelpermMultiplierColumn(); + const auto& kro = msfnTable.getOilRelpermMultiplierColumn(); + + mkrsg_[regionIdx] = NonuniformTableLinear(Ssg, krsg); + mkro_[regionIdx] = NonuniformTableLinear(Ssg, kro); + + } + } + + const TableContainer& sorwmisTables = tables->getSorwmisTables(); + if (!sorwmisTables.empty()) { + + int numRegions = sorwmisTables.size(); + + if(numRegions > 1) { + OPM_THROW(std::runtime_error, "Only single table miscibility function supported for SORWMIS"); + } + // resize the attributes of the object + sorwmis_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + const Opm::SorwmisTable& sorwmisTable = sorwmisTables.getTable(regionIdx); + + // Copy data + const auto& sw = sorwmisTable.getWaterSaturationColumn(); + const auto& sorwmis = sorwmisTable.getMiscibleResidualOilColumn(); + + sorwmis_[regionIdx] = NonuniformTableLinear(sw, sorwmis); + } + } + + const TableContainer& sgcwmisTables = tables->getSgcwmisTables(); + if (!sgcwmisTables.empty()) { + + int numRegions = sgcwmisTables.size(); + + if(numRegions > 1) { + OPM_THROW(std::runtime_error, "Only single table miscibility function supported for SGCWMIS"); + } + // resize the attributes of the object + sgcwmis_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + const Opm::SgcwmisTable& sgcwmisTable = sgcwmisTables.getTable(regionIdx); + + // Copy data + const auto& sw = sgcwmisTable.getWaterSaturationColumn(); + const auto& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn(); + + sgcwmis_[regionIdx] = NonuniformTableLinear(sw, sgcwmis); + } + } + + if (deck->hasKeyword("TLMIXPAR")) { + const auto tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(0); + std::vector mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData(); + const int numRegions = mix_params_viscosity.size(); + if (numRegions > 1) { + OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR."); + } + mix_param_viscosity_ = mix_params_viscosity[0]; + + std::vector mix_params_density = tlmixparRecord->getItem("TL_DENSITY_PARAMETER")->getSIDoubleData(); + const int numDensityItems = mix_params_density.size(); + if (numDensityItems == 0) { + mix_param_density_ = mix_param_viscosity_; + } else if (numDensityItems == 1) { + mix_param_density_ = mix_params_density[0]; + } else { + OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR."); + } + } else { + mix_param_viscosity_ = 0.0; + mix_param_density_ = 0.0; + } + + + + } } + } ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg, @@ -149,75 +300,104 @@ ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg, } ADB SolventPropsAdFromDeck::bSolvent(const ADB& pg, - const Cells& cells) const + const Cells& cells) const { - const int n = cells.size(); - assert(pg.size() == n); + return SolventPropsAdFromDeck::makeADBfromTables(pg, cells, b_); - V b(n); - V dbdp(n); - for (int i = 0; i < n; ++i) { - const double& pg_i = pg.value()[i]; - int regionIdx = cellPvtRegionIdx_[i]; - b[i] = b_[regionIdx](pg_i); - dbdp[i] = b_[regionIdx].derivative(pg_i); - } - - ADB::M dbdp_diag(dbdp.matrix().asDiagonal()); - const int num_blocks = pg.numBlocks(); - std::vector jacs(num_blocks); - for (int block = 0; block < num_blocks; ++block) { - fastSparseProduct(dbdp_diag, pg.derivative()[block], jacs[block]); - } - return ADB::function(std::move(b), std::move(jacs)); } ADB SolventPropsAdFromDeck::gasRelPermMultiplier(const ADB& solventFraction, - const Cells& cells) const + const Cells& cells) const { - const int n = cells.size(); - assert(solventFraction.value().size() == n); - V krg(n); - V dkrgdsf(n); - for (int i = 0; i < n; ++i) { - const double& solventFraction_i = solventFraction.value()[i]; - int regionIdx = 0; // TODO add mapping from cells to sat function table - krg[i] = krg_[regionIdx](solventFraction_i); - dkrgdsf[i] = krg_[regionIdx].derivative(solventFraction_i); - } + return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krg_); - ADB::M dkrgdsf_diag(dkrgdsf.matrix().asDiagonal()); - const int num_blocks = solventFraction.numBlocks(); - std::vector jacs(num_blocks); - for (int block = 0; block < num_blocks; ++block) { - fastSparseProduct(dkrgdsf_diag, solventFraction.derivative()[block], jacs[block]); - } - return ADB::function(std::move(krg), std::move(jacs)); } ADB SolventPropsAdFromDeck::solventRelPermMultiplier(const ADB& solventFraction, - const Cells& cells) const + const Cells& cells) const { + return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krs_); +} + + +ADB SolventPropsAdFromDeck::misicibleHydrocarbonWaterRelPerm(const ADB& Sn, + const Cells& cells) const +{ + return SolventPropsAdFromDeck::makeADBfromTables(Sn, cells, krn_); +} + +ADB SolventPropsAdFromDeck::miscibleSolventGasRelPermMultiplier(const ADB& Ssg, + const Cells& cells) const +{ + if (mkrsg_.size() > 0) { + return SolventPropsAdFromDeck::makeADBfromTables(Ssg, cells, mkrsg_); + } + // trivial function if not specified + return Ssg; +} + +ADB SolventPropsAdFromDeck::miscibleOilRelPermMultiplier(const ADB& So, + const Cells& cells) const +{ + if (mkro_.size() > 0) { + return SolventPropsAdFromDeck::makeADBfromTables(So, cells, mkro_); + } + // trivial function if not specified + return So; +} + +ADB SolventPropsAdFromDeck::miscibilityFunction(const ADB& solventFraction, + const Cells& cells) const +{ + + return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, misc_); +} + + +ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw, + const Cells& cells) const { + if (sgcwmis_.size()>0) { + return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sgcwmis_); + } + // return zeros if not specified + return ADB::constant(V::Zero(Sw.size())); +} + + +ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw, + const Cells& cells) const { + if (sorwmis_.size()>0) { + return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sorwmis_); + } + // return zeros if not specified + return ADB::constant(V::Zero(Sw.size())); +} + +ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD, + const Cells& cells, + const std::vector>& tables) const { const int n = cells.size(); - assert(solventFraction.value().size() == n); - V krs(n); - V dkrsdsf(n); + assert(X_AD.value().size() == n); + V x(n); + V dx(n); for (int i = 0; i < n; ++i) { - const double& solventFraction_i = solventFraction.value()[i]; + const double& X_i = X_AD.value()[i]; int regionIdx = 0; // TODO add mapping from cells to sat function table - krs[i] = krs_[regionIdx](solventFraction_i); - dkrsdsf[i] = krs_[regionIdx].derivative(solventFraction_i); + x[i] = tables[regionIdx](X_i); + dx[i] = tables[regionIdx].derivative(X_i); } - ADB::M dkrsdsf_diag(dkrsdsf.matrix().asDiagonal()); - const int num_blocks = solventFraction.numBlocks(); + ADB::M dx_diag(dx.matrix().asDiagonal()); + const int num_blocks = X_AD.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { - fastSparseProduct(dkrsdsf_diag, solventFraction.derivative()[block], jacs[block]); + fastSparseProduct(dx_diag, X_AD.derivative()[block], jacs[block]); } - return ADB::function(std::move(krs), std::move(jacs)); + return ADB::function(std::move(x), std::move(jacs)); } + + V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const { const int n = cells.size(); V density(n); @@ -228,4 +408,13 @@ V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const { return density; } +double SolventPropsAdFromDeck::mixingParameterViscosity() const { + return mix_param_viscosity_; +} + +double SolventPropsAdFromDeck::mixingParameterDensity() const { + return mix_param_density_; +} + + } //namespace OPM diff --git a/opm/autodiff/SolventPropsAdFromDeck.hpp b/opm/autodiff/SolventPropsAdFromDeck.hpp index 6cff92070..5cc4d0dab 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.hpp +++ b/opm/autodiff/SolventPropsAdFromDeck.hpp @@ -56,35 +56,95 @@ public: /// \param[in] cells Array of n cell indices to be associated with the pressure values. /// \return Array of n formation volume factor values. ADB bSolvent(const ADB& pg, - const Cells& cells) const; + const Cells& cells) const; /// Solvent viscosity. /// \param[in] pg Array of n gas pressure values. /// \param[in] cells Array of n cell indices to be associated with the pressure values. /// \return Array of n viscosity values. ADB muSolvent(const ADB& pg, - const Cells& cells) const; + const Cells& cells) const; /// Gas relPerm multipliers - /// \param[in] solventFraction Array of n solvent fraction values. + /// \param[in] gasFraction Array of n gas fraction Sg / (sg + Ss) values. /// \param[in] cells Array of n cell indices to be associated with the fraction values. /// \return Array of n gas relPerm multiplier values. ADB gasRelPermMultiplier(const ADB& solventFraction, - const Cells& cells) const; + const Cells& cells) const; /// Solvent relPerm multipliers - /// \param[in] solventFraction Array of n solvent fraction values. + /// \param[in] solventFraction Array of n solvent fraction Ss / (Sg + Ss) values. /// \param[in] cells Array of n cell indices to be associated with the fraction values. /// \return Array of n solvent relPerm multiplier values. ADB solventRelPermMultiplier(const ADB& solventFraction, - const Cells& cells) const; + const Cells& cells) const; + + /// Miscible hydrocrabon relPerm wrt water + /// \param[in] Sn Array of n total hyrdrocarbon saturation values. + /// \param[in] cells Array of n cell indices to be associated with the fraction values. + /// \return Array of n miscible hyrdrocabon wrt water relPerm values. + ADB misicibleHydrocarbonWaterRelPerm(const ADB& Sn, + const Cells& cells) const; + + /// Miscible Solvent + Gas relPerm multiplier + /// \param[in] Ssg Array of n total gas fraction (Sgas + Ssolvent) / Sn values, where + /// Sn = Sgas + Ssolvent + Soil. + /// \param[in] cells Array of n cell indices to be associated with the fraction values. + /// \return Array of n solvent gas relperm multiplier. + ADB miscibleSolventGasRelPermMultiplier(const ADB& Ssg, + const Cells& cells) const; + + /// Miscible Oil relPerm multiplier + /// \param[in] So Array of n oil fraction values. Soil / Sn values, where Sn = Sgas + Ssolvent + Soil. + /// \param[in] cells Array of n cell indices to be associated with the fraction values. + /// \return Array of n oil relperm multiplier. + ADB miscibleOilRelPermMultiplier(const ADB& So, + const Cells& cells) const; + + /// Miscible function + /// \param[in] solventFraction Array of n solvent fraction Ss / (Sg + Ss) values. + /// \param[in] cells Array of n cell indices to be associated with the fraction values. + /// \return Array of n miscibility values + ADB miscibilityFunction(const ADB& solventFraction, + const Cells& cells) const; + + /// Miscible critical gas saturation function + /// \param[in] Sw Array of n water saturation values. + /// \param[in] cells Array of n cell indices to be associated with the saturation values. + /// \return Array of n miscible critical gas saturation values + ADB miscibleCriticalGasSaturationFunction(const ADB& Sw, + const Cells& cells) const; + + /// Miscible residual oil saturation function + /// \param[in] Sw Array of n water saturation values. + /// \param[in] cells Array of n cell indices to be associated with the saturation values. + /// \return Array of n miscible residual oil saturation values + ADB miscibleResidualOilSaturationFunction(const ADB& Sw, + const Cells& cells) const; /// Solvent surface density /// \param[in] cells Array of n cell indices to be associated with the fraction values. /// \return Array of n solvent density values. V solventSurfaceDensity(const Cells& cells) const; + /// Todd-Longstaff mixing parameter for viscosity calculation + double mixingParameterViscosity() const; + + /// Todd-Longstaff mixing parameter for density calculation + double mixingParameterDensity() const; + + private: + + /// Makes ADB from table values + /// \param[in] X Array of n table lookup values. + /// \param[in] cells Array of n cell indices to be associated with the fraction values. + /// \param[in] tables Vector of tables, one for each PVT region. + /// \return Array of n solvent density values. + ADB makeADBfromTables(const ADB& X, + const Cells& cells, + const std::vector>& tables) const; + // The PVT region which is to be used for each cell std::vector cellPvtRegionIdx_; std::vector > b_; @@ -93,6 +153,14 @@ private: std::vector solvent_surface_densities_; std::vector > krg_; std::vector > krs_; + std::vector > krn_; + std::vector > mkro_; + std::vector > mkrsg_; + std::vector > misc_; + std::vector > sorwmis_; + std::vector > sgcwmis_; + double mix_param_viscosity_; + double mix_param_density_; }; } // namespace OPM diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 1a2583636..53d452e21 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -238,6 +238,8 @@ namespace Opm { computeMassFlux(const int actph , const V& transi, const ADB& kr , + const ADB& mu , + const ADB& rho , const ADB& p , const SolutionState& state ); diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index fe93b4ed0..83dead022 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -300,7 +300,10 @@ namespace Opm { } for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) { - computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state); + const std::vector& cond = phaseCondition(); + const ADB mu = fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond); + const ADB rho = fluidDensity(canph_[phaseIdx], rq_[phaseIdx].b, state.rs, state.rv); + computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], mu, rho, state.canonical_phase_pressures[canph_[phaseIdx]], state); residual_.material_balance_eq[ phaseIdx ] = pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0]) @@ -413,18 +416,18 @@ namespace Opm { BlackoilPolymerModel::computeMassFlux(const int actph , const V& transi, const ADB& kr , + const ADB& mu , + const ADB& rho , const ADB& phasePressure, const SolutionState& state) { - Base::computeMassFlux(actph, transi, kr, phasePressure, state); + Base::computeMassFlux(actph, transi, kr, mu, rho, phasePressure, state); // Polymer treatment. const int canonicalPhaseIdx = canph_[ actph ]; if (canonicalPhaseIdx == Water) { if (has_polymer_) { - const std::vector& cond = phaseCondition(); const ADB tr_mult = transMult(state.pressure); - const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv, cond); const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); const ADB mc = computeMc(state); const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr);