From dcc6eb2e0cd37873e539502809634ceb3d82718f Mon Sep 17 00:00:00 2001
From: Paul Egberts <paul.egberts@tno.nl>
Date: Wed, 18 May 2022 16:01:10 +0200
Subject: [PATCH] modified density calculation to include oil vapor and
 dissolved gas

---
 .../fluidsystems/blackoilpvt/GasPvtThermal.hpp | 18 +++++++++++++++---
 .../fluidsystems/blackoilpvt/OilPvtThermal.hpp | 16 ++++++++++++++--
 2 files changed, 29 insertions(+), 5 deletions(-)

diff --git a/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp b/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp
index ae7c64a57..479d4cea9 100644
--- a/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp
+++ b/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp
@@ -160,6 +160,13 @@ public:
                 gasJTRefPres_[regionIdx] =  record.P0;
                 gasJTC_[regionIdx] = record.C1;
             }
+
+            const auto& densityTable = eclState.getTableManager().getDensityTable();
+
+            assert(densityTable.size() == numRegions);
+            for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
+                 rhoRefO_[regionIdx] = densityTable[regionIdx].oil;
+            }
         }
 
         if (enableInternalEnergy_) {
@@ -214,6 +221,7 @@ public:
         gasdentCT2_.resize(numRegions);
         gasJTRefPres_.resize(numRegions);
         gasJTC_.resize(numRegions);
+        rhoRefO_.resize(numRegions);
     }
 
     /*!
@@ -269,7 +277,7 @@ public:
             Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rv);
             const Scalar hVap = 480.6e3; // [J / kg]
             Evaluation Cp = (internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true) - hVap)/temperature;
-            Evaluation density = invB * gasReferenceDensity(regionIdx);
+            Evaluation density = invB * (gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
 
             Evaluation enthalpyPres;
             if  (JTC != 0) {
@@ -287,8 +295,10 @@ public:
                 Evaluation enthalpyPresPrev = 0;
                 for (size_t i = 0; i < N; ++i) {
                     Evaluation Pnew = Pref + i * deltaP;
-                    Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, Rv) * gasReferenceDensity(regionIdx);
-                    Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;  
+                    Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, Rv) *
+                                     (gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
+                    // see e.g.https://en.wikipedia.org/wiki/Joule-Thomson_effect for a derivation of the Joule-Thomson coeff.
+                    Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
                     Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
                     enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres; 
                     enthalpyPresPrev = enthalpyPres;
@@ -545,6 +555,8 @@ private:
     std::vector<Scalar> gasJTRefPres_;
     std::vector<Scalar> gasJTC_;
 
+    std::vector<Scalar> rhoRefO_;
+
     // piecewise linear curve representing the internal energy of gas
     std::vector<TabulatedOneDFunction> internalEnergyCurves_;
 
diff --git a/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.hpp b/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.hpp
index 2dfc01a36..f98f25743 100644
--- a/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.hpp
+++ b/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.hpp
@@ -185,6 +185,13 @@ public:
                 oilJTRefPres_[regionIdx] =  record.P0;
                 oilJTC_[regionIdx] = record.C1;
             }
+
+            const auto& densityTable = eclState.getTableManager().getDensityTable();
+
+            assert(densityTable.size() == numRegions);
+            for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
+                 rhoRefG_[regionIdx] = densityTable[regionIdx].gas;
+            }
         }
 
         if (enableInternalEnergy_) {
@@ -235,6 +242,7 @@ public:
         oildentCT2_.resize(numRegions);
         oilJTRefPres_.resize(numRegions);
         oilJTC_.resize(numRegions);
+        rhoRefG_.resize(numRegions);
     }
 
     /*!
@@ -289,7 +297,7 @@ public:
 
             Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rs);
             Evaluation Cp = internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true)/temperature;
-            Evaluation density = invB * oilReferenceDensity(regionIdx);
+            Evaluation density = invB * (oilReferenceDensity(regionIdx) + Rs * rhoRefG_[regionIdx]);
 
             Evaluation enthalpyPres;
             if  (JTC != 0) {
@@ -307,7 +315,9 @@ public:
                 Evaluation enthalpyPresPrev = 0;
                 for (size_t i = 0; i < N; ++i) {
                     Evaluation Pnew = Pref + i * deltaP;
-                    Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, Rs) * oilReferenceDensity(regionIdx);
+                    Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, Rs) *
+                                     (oilReferenceDensity(regionIdx) + Rs * rhoRefG_[regionIdx]) ;
+                    // see e.g.https://en.wikipedia.org/wiki/Joule-Thomson_effect for a derivation of the Joule-Thomson coeff.
                     Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;  
                     Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
                     enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres; 
@@ -564,6 +574,8 @@ private:
     std::vector<Scalar> oilJTRefPres_;
     std::vector<Scalar> oilJTC_;
 
+    std::vector<Scalar> rhoRefG_;
+
     // piecewise linear curve representing the internal energy of oil
     std::vector<TabulatedOneDFunction> internalEnergyCurves_;