Skip to content

water viscosity calculation for brine + thermal #4656

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
plgbrts wants to merge 1 commit into
base: master
Choose a base branch
from
Open

water viscosity calculation for brine + thermal #4656

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 30 additions & 12 deletions opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
enableJouleThomson_ = tables.WatJT().size() > 0;
enableThermalViscosity_ = tables.hasTables("WATVISCT");
enableInternalEnergy_ = tables.hasTables("SPECHEAT");
enableBrineViscosity_ = eclState.runspec().phases().active(Phase::BRINE);

unsigned regions = isothermalPvt_->numRegions();
setNumRegions(regions);
Expand Down Expand Up @@ -115,13 +116,30 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
const auto& watvisctTables = tables.getWatvisctTables();
const auto& viscrefTables = tables.getViscrefTable();

const auto& pvtwTables = tables.getPvtwTable();

if (pvtwTables.size() != regions) {
OPM_THROW(std::runtime_error,
fmt::format("Table sizes mismatch. PVTW: {}, numRegions: {}\n",
pvtwTables.size(), regions));
if (enableBrineViscosity_) {
const auto& pvtwSaltTables = tables.getPvtwSaltTables();
if (pvtwSaltTables.size() != regions) {
OPM_THROW(std::runtime_error,
fmt::format("Table sizes mismatch. PVTWSALT: {}, numRegions: {}\n",
pvtwSaltTables.size(), regions));
}
for (unsigned regionIdx = 0; regionIdx < regions; ++ regionIdx) {
referenceSaltConcentration_[regionIdx] = pvtwSaltTables[regionIdx].getReferenceSaltConcentrationValue();
}
Copy link
Member

@totto82 totto82 Nov 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You will need to set the viscosity and viscosibility also for PVTWSALT

Copy link
Contributor Author

@plgbrts plgbrts Jan 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This indeed solved the assertion failure I got after a rebase. Thank for pointing out this.

}
else {
const auto& pvtwTables = tables.getPvtwTable();
if (pvtwTables.size() != regions) {
OPM_THROW(std::runtime_error,
fmt::format("Table sizes mismatch. PVTW: {}, numRegions: {}\n",
pvtwTables.size(), regions));
}
for (unsigned regionIdx = 0; regionIdx < regions; ++ regionIdx) {
pvtwViscosity_[regionIdx] = pvtwTables[regionIdx].viscosity;
pvtwViscosibility_[regionIdx] = pvtwTables[regionIdx].viscosibility;
}
}

if (watvisctTables.size() != regions) {
OPM_THROW(std::runtime_error,
fmt::format("Table sizes mismatch. WATVISCT: {}, numRegions: {}\n",
Expand All @@ -140,11 +158,6 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)

viscrefPress_[regionIdx] = viscrefTables[regionIdx].reference_pressure;
}

for (unsigned regionIdx = 0; regionIdx < regions; ++ regionIdx) {
pvtwViscosity_[regionIdx] = pvtwTables[regionIdx].viscosity;
pvtwViscosibility_[regionIdx] = pvtwTables[regionIdx].viscosibility;
}
}

if (enableInternalEnergy_) {
Expand Down Expand Up @@ -189,6 +202,7 @@ setNumRegions(std::size_t numRegions)
pvtwCompressibility_.resize(numRegions);
pvtwViscosity_.resize(numRegions);
pvtwViscosibility_.resize(numRegions);
referenceSaltConcentration_.resize(numRegions);
viscrefPress_.resize(numRegions);
watvisctCurves_.resize(numRegions);
watdentRefTemp_.resize(numRegions);
Expand Down Expand Up @@ -222,12 +236,14 @@ operator==(const WaterPvtThermal<Scalar, enableBrine>& data) const
this->pvtwCompressibility() == data.pvtwCompressibility() &&
this->pvtwViscosity() == data.pvtwViscosity() &&
this->pvtwViscosibility() == data.pvtwViscosibility() &&
this->referenceSaltConcentration() == data.referenceSaltConcentration() &&
this->watvisctCurves() == data.watvisctCurves() &&
this->internalEnergyCurves() == data.internalEnergyCurves() &&
this->enableThermalDensity() == data.enableThermalDensity() &&
this->enableJouleThomson() == data.enableJouleThomson() &&
this->enableThermalViscosity() == data.enableThermalViscosity() &&
this->enableInternalEnergy() == data.enableInternalEnergy();
this->enableInternalEnergy() == data.enableInternalEnergy() &&
this->enableBrineViscosity() == data.enableBrineViscosity();
}

template<class Scalar, bool enableBrine>
Expand All @@ -252,12 +268,14 @@ operator=(const WaterPvtThermal<Scalar, enableBrine>& data)
pvtwCompressibility_ = data.pvtwCompressibility_;
pvtwViscosity_ = data.pvtwViscosity_;
pvtwViscosibility_ = data.pvtwViscosibility_;
referenceSaltConcentration_ = data.referenceSaltConcentration_;
watvisctCurves_ = data.watvisctCurves_;
internalEnergyCurves_ = data.internalEnergyCurves_;
enableThermalDensity_ = data.enableThermalDensity_;
enableJouleThomson_ = data.enableJouleThomson_;
enableThermalViscosity_ = data.enableThermalViscosity_;
enableInternalEnergy_ = data.enableInternalEnergy_;
enableBrineViscosity_ = data.enableBrineViscosity_;

return *this;
}
Expand Down
30 changes: 22 additions & 8 deletions opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,12 +68,14 @@ class WaterPvtThermal
const std::vector<Scalar>& pvtwCompressibility,
const std::vector<Scalar>& pvtwViscosity,
const std::vector<Scalar>& pvtwViscosibility,
const std::vector<Scalar>& referenceSaltConcentration,
const std::vector<TabulatedOneDFunction>& watvisctCurves,
const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
bool enableThermalDensity,
bool enableJouleThomson,
bool enableThermalViscosity,
bool enableInternalEnergy)
bool enableInternalEnergy,
bool enableBrineViscosity)
: isothermalPvt_(isothermalPvt)
, viscrefPress_(viscrefPress)
, watdentRefTemp_(watdentRefTemp)
Expand All @@ -86,12 +88,14 @@ class WaterPvtThermal
, pvtwCompressibility_(pvtwCompressibility)
, pvtwViscosity_(pvtwViscosity)
, pvtwViscosibility_(pvtwViscosibility)
, referenceSaltConcentration_(referenceSaltConcentration)
, watvisctCurves_(watvisctCurves)
, internalEnergyCurves_(internalEnergyCurves)
, enableThermalDensity_(enableThermalDensity)
, enableJouleThomson_(enableJouleThomson)
, enableThermalViscosity_(enableThermalViscosity)
, enableInternalEnergy_(enableInternalEnergy)
, enableBrineViscosity_(enableBrineViscosity)
{ }

WaterPvtThermal(const WaterPvtThermal& data)
Expand Down Expand Up @@ -227,14 +231,16 @@ class WaterPvtThermal
if (!enableThermalViscosity()) {
return isothermalMu;
}

Scalar x = -pvtwViscosibility_[regionIdx] * (viscrefPress_[regionIdx] -
pvtwRefPress_[regionIdx]);
Scalar muRef = pvtwViscosity_[regionIdx] / (1.0 + x + 0.5 * x * x);

// compute the viscosity deviation due to temperature

const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
return isothermalMu * muWatvisct / muRef;
if (enableBrineViscosity()) {
auto muRef = isothermalPvt_->viscosity(regionIdx, temperature, Evaluation(viscrefPress_[regionIdx]), Rsw, Evaluation(referenceSaltConcentration_[regionIdx]));
return isothermalMu * muWatvisct / muRef;
} else {
Scalar x = -pvtwViscosibility_[regionIdx] * (viscrefPress_[regionIdx] - pvtwRefPress_[regionIdx]);
Scalar muRef = pvtwViscosity_[regionIdx] / (1.0 + x + 0.5 * x * x);
return isothermalMu * muWatvisct / muRef;
}
}

/*!
Expand Down Expand Up @@ -404,6 +410,9 @@ class WaterPvtThermal
const std::vector<Scalar>& pvtwViscosibility() const
{ return pvtwViscosibility_; }

const std::vector<Scalar>& referenceSaltConcentration() const
{ return referenceSaltConcentration_; }

const std::vector<TabulatedOneDFunction>& watvisctCurves() const
{ return watvisctCurves_; }

Expand All @@ -413,6 +422,9 @@ class WaterPvtThermal
bool enableInternalEnergy() const
{ return enableInternalEnergy_; }

bool enableBrineViscosity() const
{ return enableBrineViscosity_; }

const std::vector<Scalar>& watJTRefPres() const
{ return watJTRefPres_; }

Expand Down Expand Up @@ -443,6 +455,7 @@ class WaterPvtThermal
std::vector<Scalar> pvtwCompressibility_{};
std::vector<Scalar> pvtwViscosity_{};
std::vector<Scalar> pvtwViscosibility_{};
std::vector<Scalar> referenceSaltConcentration_{};

std::vector<TabulatedOneDFunction> watvisctCurves_{};

Expand All @@ -454,6 +467,7 @@ class WaterPvtThermal
bool enableJouleThomson_{false};
bool enableThermalViscosity_{false};
bool enableInternalEnergy_{false};
bool enableBrineViscosity_{false};
};

} // namespace Opm
Expand Down