From 2d166c7ce0d17cc71539f4ad2af754e9c1cd703e Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Fri, 27 Apr 2018 14:27:27 +0200
Subject: [PATCH] [fluidstates] make it possible to use fluidstates with
 several temperatures

---
 dumux/material/fluidstates/compositional.hh  | 15 +++--
 dumux/material/fluidstates/immiscible.hh     | 20 +++++--
 dumux/material/fluidstates/nonequilibrium.hh | 61 ++++++++++++++++++++
 3 files changed, 85 insertions(+), 11 deletions(-)

diff --git a/dumux/material/fluidstates/compositional.hh b/dumux/material/fluidstates/compositional.hh
index 2e065d34d1..f1eed0374f 100644
--- a/dumux/material/fluidstates/compositional.hh
+++ b/dumux/material/fluidstates/compositional.hh
@@ -208,7 +208,7 @@ public:
      * \brief The absolute temperature\f$T_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[K]}\f$
      */
     Scalar temperature(int phaseIdx) const
-    { return temperature_; }
+    { return temperature_[phaseIdx]; }
 
     /*!
      * \brief The pressure \f$p_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[Pa]}\f$
@@ -257,7 +257,7 @@ public:
      * \brief The temperature within the domain \f$\mathrm{[K]}\f$
      */
     Scalar temperature() const
-    { return temperature_; }
+    { return temperature_[0]; }
 
     /*!
      * \brief The fugacity of a component  \f$\mathrm{[Pa]}\f$
@@ -288,6 +288,7 @@ public:
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             averageMolarMass_[phaseIdx] = 0;
             sumMoleFractions_[phaseIdx] = 0;
+            temperature_[phaseIdx] = fs.temperature();
             for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
                 moleFraction_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
                 fugacityCoefficient_[phaseIdx][compIdx] = fs.fugacityCoefficient(phaseIdx, compIdx);
@@ -300,7 +301,6 @@ public:
             enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
         }
-        temperature_ = fs.temperature(0);
         wPhaseIdx_ = fs.wettingPhase();
     }
 
@@ -308,7 +308,10 @@ public:
      * \brief Set the temperature \f$\mathrm{[K]}\f$ of all phases.
      */
     void setTemperature(Scalar value)
-    { temperature_ = value; }
+    {
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+            temperature_[phaseIdx] = value;
+    }
 
     /*!
      * \brief Set the temperature \f$\mathrm{[K]}\f$ of a specific phase.
@@ -316,7 +319,7 @@ public:
      */
     void setTemperature(const int phaseIdx, const Scalar value)
     {
-        DUNE_THROW(Dune::NotImplemented, "This is a fluidstate for equilibrium, temperature in all phases is assumed to be equal.");
+       temperature_[phaseIdx] = value;
     }
 
     /*!
@@ -481,7 +484,7 @@ protected:
     Scalar density_[numPhases];
     Scalar enthalpy_[numPhases];
     Scalar viscosity_[numPhases];
-    Scalar temperature_;
+    Scalar temperature_[numPhases];
 
     // For porous medium flow models, here we ... the index of the wetting
     // phase (needed for vapor pressure evaluation if kelvin equation is used)
diff --git a/dumux/material/fluidstates/immiscible.hh b/dumux/material/fluidstates/immiscible.hh
index 8688eb646f..4496472eb1 100644
--- a/dumux/material/fluidstates/immiscible.hh
+++ b/dumux/material/fluidstates/immiscible.hh
@@ -211,7 +211,7 @@ public:
      * \brief The absolute temperature\f$T_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[K]}\f$
      */
     Scalar temperature(int phaseIdx) const
-    { return temperature_; }
+    { return temperature_[phaseIdx]; }
 
     /*!
      * \brief The pressure \f$p_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[Pa]}\f$
@@ -250,7 +250,7 @@ public:
      * \brief The temperature within the domain \f$\mathrm{[K]}\f$
      */
     Scalar temperature() const
-    { return temperature_; }
+    { return temperature_[0]; }
 
     /*!
      * \brief The fugacity of a component  \f$\mathrm{[Pa]}\f$
@@ -285,15 +285,25 @@ public:
             density_[phaseIdx] = fs.density(phaseIdx);
             enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
+            temperature_[phaseIdx] = fs.temperature(0);
         }
-        temperature_ = fs.temperature(0);
+
     }
 
+    /*!
+     * \brief Set the temperature \f$\mathrm{[K]}\f$ of a fluid phase
+     */
+    void setTemperature(int phaseIdx, Scalar value)
+    { temperature_[phaseIdx] = value; }
+
     /*!
      * \brief Set the temperature \f$\mathrm{[K]}\f$ of a fluid phase
      */
     void setTemperature(Scalar value)
-    { temperature_ = value; }
+    {
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+            temperature_[phaseIdx] = value;
+    }
 
     /*!
      * \brief Set the fluid pressure of a phase \f$\mathrm{[Pa]}\f$
@@ -357,7 +367,7 @@ protected:
     Scalar density_[numPhases];
     Scalar enthalpy_[numPhases];
     Scalar viscosity_[numPhases];
-    Scalar temperature_;
+    Scalar temperature_[numPhases];
 };
 
 } // end namespace Dumux
diff --git a/dumux/material/fluidstates/nonequilibrium.hh b/dumux/material/fluidstates/nonequilibrium.hh
index d90281815d..0fa62dea73 100644
--- a/dumux/material/fluidstates/nonequilibrium.hh
+++ b/dumux/material/fluidstates/nonequilibrium.hh
@@ -68,6 +68,13 @@ public:
      * Generic access to fluid properties (No assumptions
      * on thermodynamic equilibrium required)
      *****************************************************/
+    /*!
+     * \brief Returns the index of the wetting phase in the
+     *        fluid-solid configuration (for porous medium systems).
+     *
+     * \param phaseIdx the index of the phase
+     */
+    int wettingPhase() const { return wPhaseIdx_; }
     /*!
      * \brief Returns the saturation \f$S_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[-]}\f$.
      *
@@ -221,6 +228,15 @@ public:
     Scalar pressure(int phaseIdx) const
     { return pressure_[phaseIdx]; }
 
+    /*!
+     * \brief The partial pressure of a component in a phase \f$\mathrm{[Pa]}\f$
+     */
+    Scalar partialPressure(int phaseIdx, int compIdx) const
+    {
+        assert(FluidSystem::isGas(phaseIdx));
+        return moleFraction(phaseIdx, compIdx) * pressure(phaseIdx);
+    }
+
     /*!
      * \brief The specific enthalpy \f$h_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[J/kg]}\f$
      */
@@ -311,6 +327,41 @@ public:
         }
     }
 
+    /*!
+     * \brief Set the mass fraction of a component in a phase \f$\mathrm{[-]}\f$
+     *        and update the average molar mass \f$\mathrm{[kg/mol]}\f$ according
+     *        to the current composition of the phase
+     */
+    void setMassFraction(int phaseIdx, int compIdx, Scalar value)
+    {
+        Valgrind::CheckDefined(value);
+        Valgrind::SetDefined(sumMoleFractions_[phaseIdx]);
+        Valgrind::SetDefined(averageMolarMass_[phaseIdx]);
+        Valgrind::SetDefined(moleFraction_[phaseIdx][compIdx]);
+
+        if (numComponents != 2)
+            DUNE_THROW(Dune::NotImplemented, "This currently only works for 2 components.");
+        else
+        {
+            // calculate average molar mass of the gas phase
+            Scalar M1 = FluidSystem::molarMass(compIdx);
+            Scalar M2 = FluidSystem::molarMass(1-compIdx);
+            Scalar X2 = 1.0-value;
+            Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
+
+            moleFraction_[phaseIdx][compIdx] = value * avgMolarMass / M1;
+            moleFraction_[phaseIdx][1-compIdx] = 1.0-moleFraction_[phaseIdx][compIdx];
+
+            // re-calculate the mean molar mass
+            sumMoleFractions_[phaseIdx] = 0.0;
+            averageMolarMass_[phaseIdx] = 0.0;
+            for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
+                sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compJIdx];
+                averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compJIdx]*FluidSystem::molarMass(compJIdx);
+            }
+        }
+    }
+
     /*!
      * \brief Set the fugacity of a component in a phase \f$\mathrm{[-]}\f$
      */
@@ -335,6 +386,12 @@ public:
     void setViscosity(int phaseIdx, Scalar value)
     { viscosity_[phaseIdx] = value; }
 
+    /*!
+     * \brief Set the index of the wetting phase
+     */
+    void setWettingPhase(int phaseIdx)
+    { wPhaseIdx_ = phaseIdx; }
+
     /*!
      * \brief Retrieve all parameters from an arbitrary fluid
      *        state.
@@ -402,6 +459,10 @@ protected:
     Scalar viscosity_[numPhases];
     Scalar temperature_[numPhases]; //
     Scalar temperatureEquil_;
+
+    // For porous medium flow models, here we ... the index of the wetting
+    // phase (needed for vapor pressure evaluation if kelvin equation is used)
+    int wPhaseIdx_{0};
 };
 
 } // end namespace Dumux
-- 
GitLab