guard against d = 0

This commit is contained in:
Tor Harald Sandve 2022-03-30 08:25:35 +02:00
parent 2e91d2f353
commit 5b53fcd8a6
3 changed files with 13 additions and 10 deletions

View File

@ -758,10 +758,10 @@ computeSegmentFluidProperties(const EvalWell& temperature,
const EvalWell d = 1.0 - rs * rv;
if (rs != 0.0) { // rs > 0.0?
if (rs > 0.0 && d > 0.0) {
mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d;
}
if (rv != 0.0) { // rv > 0.0?
if (rv > 0.0 && d > 0.0) {
mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d;
}
}

View File

@ -950,13 +950,14 @@ computeConnectionDensities(const std::vector<double>& perfComponentRates,
if (!rvmax_perf.empty() && mix[gaspos] > 1e-12) {
rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]);
}
if (rs != 0.0) {
double d = 1.0 - rs*rv;
if (rs > 0.0 && d > 0.0) {
// Subtract gas in oil from gas mixture
x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv);
x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/d;
}
if (rv != 0.0) {
if (rv > 0.0 && d > 0.0) {
// Subtract oil in gas from oil mixture
x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);;
x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/d;
}
}
double volrat = 0.0;

View File

@ -356,10 +356,14 @@ namespace Opm
const double d = 1.0 - getValue(rv) * getValue(rs);
// vaporized oil into gas
// rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
if (d > 0.0) {
perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
}
// dissolved of gas in oil
// rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
if (d > 0.0) {
perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
}
}
}
}
@ -580,8 +584,6 @@ namespace Opm
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// q_os = q_or * b_o + rv * q_gr * b_g
// q_gs = q_gr * g_g + rs * q_or * b_o
// d = 1.0 - rs * rv
const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
if(FluidSystem::gasPhaseIdx == phaseIdx)
cq_r_thermal = (cq_s[gasCompIdx] - this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );