From 29f474637b7c1615e9cc32d3144f206e138cc936 Mon Sep 17 00:00:00 2001
From: Timo Koch <timokoch@uio.no>
Date: Mon, 24 Jul 2023 19:48:50 +0200
Subject: [PATCH] [bugfix][energy] Fix energy balance by addition potential
energy term
---
dumux/porousmediumflow/2p1c/localresidual.hh | 2 +-
dumux/porousmediumflow/3p3c/localresidual.hh | 2 +-
.../3pwateroil/localresidual.hh | 2 +-
.../compositional/localresidual.hh | 2 +-
.../immiscible/localresidual.hh | 2 +-
.../nonequilibrium/thermal/localresidual.hh | 27 +++++-
.../nonisothermal/localresidual.hh | 58 +++++++----
.../localresidual_incompressible.hh | 96 -------------------
dumux/porousmediumflow/nonisothermal/model.hh | 9 +-
.../richards/localresidual.hh | 4 +-
10 files changed, 78 insertions(+), 126 deletions(-)
delete mode 100644 dumux/porousmediumflow/nonisothermal/localresidual_incompressible.hh
diff --git a/dumux/porousmediumflow/2p1c/localresidual.hh b/dumux/porousmediumflow/2p1c/localresidual.hh
index 61ecea0b3f..fe41386d16 100644
--- a/dumux/porousmediumflow/2p1c/localresidual.hh
+++ b/dumux/porousmediumflow/2p1c/localresidual.hh
@@ -59,7 +59,7 @@ public:
volVars.porosity()
* volVars.saturation(phaseIdx) * volVars.density(phaseIdx);
- EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+ EnergyLocalResidual::fluidPhaseStorage(storage, problem, scv, volVars, phaseIdx);
}
// The energy storage in the solid matrix
diff --git a/dumux/porousmediumflow/3p3c/localresidual.hh b/dumux/porousmediumflow/3p3c/localresidual.hh
index 61c6d75aa8..506c06f03f 100644
--- a/dumux/porousmediumflow/3p3c/localresidual.hh
+++ b/dumux/porousmediumflow/3p3c/localresidual.hh
@@ -96,7 +96,7 @@ public:
}
// The energy storage in the fluid phase with index phaseIdx
- EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+ EnergyLocalResidual::fluidPhaseStorage(storage, problem, scv, volVars, phaseIdx);
}
// The energy storage in the solid matrix
diff --git a/dumux/porousmediumflow/3pwateroil/localresidual.hh b/dumux/porousmediumflow/3pwateroil/localresidual.hh
index 1ad1d9679e..36467645a0 100644
--- a/dumux/porousmediumflow/3pwateroil/localresidual.hh
+++ b/dumux/porousmediumflow/3pwateroil/localresidual.hh
@@ -101,7 +101,7 @@ public:
}
//! The energy storage in the fluid phase with index phaseIdx
- EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+ EnergyLocalResidual::fluidPhaseStorage(storage, problem, scv, volVars, phaseIdx);
}
//! The energy storage in the solid matrix
diff --git a/dumux/porousmediumflow/compositional/localresidual.hh b/dumux/porousmediumflow/compositional/localresidual.hh
index 3fb810b373..307cf7f7da 100644
--- a/dumux/porousmediumflow/compositional/localresidual.hh
+++ b/dumux/porousmediumflow/compositional/localresidual.hh
@@ -106,7 +106,7 @@ public:
* volVars.saturation(phaseIdx);
//! The energy storage in the fluid phase with index phaseIdx
- EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+ EnergyLocalResidual::fluidPhaseStorage(storage, problem, scv, volVars, phaseIdx);
}
//! The energy storage in the solid matrix
diff --git a/dumux/porousmediumflow/immiscible/localresidual.hh b/dumux/porousmediumflow/immiscible/localresidual.hh
index ca05a02c90..6589f5ba8f 100644
--- a/dumux/porousmediumflow/immiscible/localresidual.hh
+++ b/dumux/porousmediumflow/immiscible/localresidual.hh
@@ -74,7 +74,7 @@ public:
* volVars.saturation(phaseIdx);
//! The energy storage in the fluid phase with index phaseIdx
- EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+ EnergyLocalResidual::fluidPhaseStorage(storage, problem, scv, volVars, phaseIdx);
}
//! The energy storage in the solid matrix
diff --git a/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh b/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
index 75c54178af..4fda272afa 100644
--- a/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
+++ b/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
@@ -37,6 +37,7 @@ class EnergyLocalResidualNonEquilibrium<TypeTag, 1/*numEnergyEqFluid*/>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -59,17 +60,27 @@ class EnergyLocalResidualNonEquilibrium<TypeTag, 1/*numEnergyEqFluid*/>
static constexpr auto numComponents = ModelTraits::numFluidComponents();
public:
+ static void fluidPhaseStorage(NumEqVector& storage,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars,
+ int phaseIdx)
+ {
+ static_assert("Deprecated interface that has been removed!");
+ }
+
//! The energy storage in the fluid phase with index phaseIdx
static void fluidPhaseStorage(NumEqVector& storage,
+ const Problem&,
const SubControlVolume& scv,
const VolumeVariables& volVars,
int phaseIdx)
{
- //in case we have one energy equation for more than one fluid phase, add up parts on the one energy equation
+ // in case we have one energy equation for more than one fluid phase,
+ // add up parts on the one energy equation
storage[energyEq0Idx] += volVars.porosity()
* volVars.density(phaseIdx)
- * volVars.internalEnergy(phaseIdx)
- * volVars.saturation(phaseIdx);
+ * volVars.saturation(phaseIdx)
+ * (volVars.internalEnergy(phaseIdx) - gravityPotential);
}
@@ -213,6 +224,7 @@ class EnergyLocalResidualNonEquilibrium<TypeTag, 2/*numEnergyEqFluid*/>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -242,9 +254,17 @@ class EnergyLocalResidualNonEquilibrium<TypeTag, 2/*numEnergyEqFluid*/>
static constexpr bool enableChemicalNonEquilibrium = ModelTraits::enableChemicalNonEquilibrium();
public:
+ static void fluidPhaseStorage(NumEqVector& storage,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars,
+ int phaseIdx)
+ {
+ static_assert("Deprecated interface that has been removed!");
+ }
//! The energy storage in the fluid phase with index phaseIdx
static void fluidPhaseStorage(NumEqVector& storage,
+ const Problem&,
const SubControlVolume& scv,
const VolumeVariables& volVars,
int phaseIdx)
@@ -253,6 +273,7 @@ public:
* volVars.density(phaseIdx)
* volVars.internalEnergy(phaseIdx)
* volVars.saturation(phaseIdx);
+
}
//! The advective phase energy fluxes
diff --git a/dumux/porousmediumflow/nonisothermal/localresidual.hh b/dumux/porousmediumflow/nonisothermal/localresidual.hh
index f203eb31f4..0d5708405d 100644
--- a/dumux/porousmediumflow/nonisothermal/localresidual.hh
+++ b/dumux/porousmediumflow/nonisothermal/localresidual.hh
@@ -8,7 +8,7 @@
* \file
* \ingroup NIModel
* \brief Element-wise calculation of the local residual for non-isothermal
- * fully implicit models.
+ * fully implicit models. See NIModel for the detailed description.
*/
#ifndef DUMUX_ENERGY_LOCAL_RESIDUAL_HH
@@ -28,28 +28,33 @@ using EnergyLocalResidual = EnergyLocalResidualImplementation<TypeTag, GetPropTy
/*!
* \ingroup NIModel
- * \brief Element-wise calculation of the energy residual for non-isothermal problems.
+ * \brief Element-wise calculation of the energy residual for isothermal problems.
*/
template<class TypeTag>
class EnergyLocalResidualImplementation<TypeTag, false>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
public:
+ static void fluidPhaseStorage(NumEqVector& storage,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars,
+ int phaseIdx)
+ {
+ static_assert("Deprecated interface that has been removed!");
+ }
+
/*!
* \brief The energy storage in the fluid phase with index phaseIdx.
- *
- * \param storage The mass of the component within the sub-control volume
- * \param scv The sub-control volume
- * \param volVars The volume variables
- * \param phaseIdx The phase index
*/
static void fluidPhaseStorage(NumEqVector& storage,
+ const Problem& problem,
const SubControlVolume& scv,
const VolumeVariables& volVars,
int phaseIdx)
@@ -102,13 +107,14 @@ public:
/*!
* \ingroup NIModel
- * \brief TODO docme!
+ * \brief Element-wise calculation of the energy residual for non-isothermal problems.
*/
template<class TypeTag>
class EnergyLocalResidualImplementation<TypeTag, true>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
@@ -125,22 +131,31 @@ class EnergyLocalResidualImplementation<TypeTag, true>
public:
/*!
- * \brief The energy storage in the fluid phase with index phaseIdx
- *
- * \param storage The mass of the component within the sub-control volume
- * \param scv The sub-control volume
- * \param volVars The volume variables
- * \param phaseIdx The phase index
+ * \brief The energy storage in the fluid phase with index phaseIdx.
*/
static void fluidPhaseStorage(NumEqVector& storage,
+ const Problem& problem,
const SubControlVolume& scv,
const VolumeVariables& volVars,
int phaseIdx)
{
+ // this implementation of the potential energy contribution is only correct
+ // if gravity vector is constant in space and time
+ const auto& x = scv.dofPosition();
+ const auto gravityPotential = x*problem.spatialParams().gravity(x);
+
storage[energyEqIdx] += volVars.porosity()
* volVars.density(phaseIdx)
- * volVars.internalEnergy(phaseIdx)
- * volVars.saturation(phaseIdx);
+ * volVars.saturation(phaseIdx)
+ * (volVars.internalEnergy(phaseIdx) - gravityPotential);
+ }
+
+ static void fluidPhaseStorage(NumEqVector& storage,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars,
+ int phaseIdx)
+ {
+ static_assert("Deprecated interface that has been removed!");
}
/*!
@@ -171,8 +186,15 @@ public:
FluxVariables& fluxVars,
int phaseIdx)
{
- auto upwindTerm = [phaseIdx](const auto& volVars)
- { return volVars.density(phaseIdx)*volVars.mobility(phaseIdx)*volVars.enthalpy(phaseIdx); };
+ // this implementation of the potential energy contribution is only correct
+ // if gravity vector is constant in space and time
+ const auto& x = fluxVars.scvFace().ipGlobal();
+ const auto gravityPotential = x*fluxVars.problem().spatialParams().gravity(x);
+
+ auto upwindTerm = [=](const auto& volVars){
+ return volVars.density(phaseIdx)*volVars.mobility(phaseIdx)
+ * (volVars.enthalpy(phaseIdx) - gravityPotential);
+ };
flux[energyEqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
}
diff --git a/dumux/porousmediumflow/nonisothermal/localresidual_incompressible.hh b/dumux/porousmediumflow/nonisothermal/localresidual_incompressible.hh
deleted file mode 100644
index 850d4b474c..0000000000
--- a/dumux/porousmediumflow/nonisothermal/localresidual_incompressible.hh
+++ /dev/null
@@ -1,96 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup NIModel
- * \brief Element-wise calculation of the local residual for non-isothermal
- * fully implicit models.
- */
-
-#ifndef DUMUX_ENERGY_LOCAL_RESIDUAL_INCOMPRESSIBLE_HH
-#define DUMUX_ENERGY_LOCAL_RESIDUAL_INCOMPRESSIBLE_HH
-
-#include <dumux/common/properties.hh>
-#include <dumux/common/numeqvector.hh>
-#include "localresidual.hh"
-
-namespace Dumux {
-
-// forward declaration
-template<class TypeTag, bool enableEneryBalance>
-class EnergyLocalResidualIncompressibleImplementation;
-
-template<class TypeTag>
-using EnergyLocalResidualIncompressible = EnergyLocalResidualIncompressibleImplementation<TypeTag, GetPropType<TypeTag, Properties::ModelTraits>::enableEnergyBalance()>;
-
-/*!
- * \ingroup NIModel
- * \brief Element-wise calculation of the energy residual for non-isothermal problems.
- */
-template<class TypeTag>
-class EnergyLocalResidualIncompressibleImplementation<TypeTag, false>: public EnergyLocalResidualImplementation<TypeTag, false>
-{}; //energy balance not enabled
-
-/*!
- * \ingroup NIModel
- * \brief TODO docme!
- */
-template<class TypeTag>
-class EnergyLocalResidualIncompressibleImplementation<TypeTag, true>: public EnergyLocalResidualImplementation<TypeTag, true>
-{
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
- using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
- using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
- using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Element = typename GridView::template Codim<0>::Entity;
- using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
- using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
- using Indices = typename ModelTraits::Indices;
-
- static constexpr int numPhases = ModelTraits::numFluidPhases();
- enum { energyEqIdx = Indices::energyEqIdx };
-
-public:
- /*!
- * \brief The advective phase energy fluxes for incompressible flow.
- *
- * Using specific internal energy $u$ instead of specific enthalpy \f$h\f$ for incompressible flow in convective flux
- * to account for otherwise neglected pressure work term (\f$\nabla p \cdot v\f$).
- *
- * Compressible formulation in EnergyLocalResidual (neglecting pressure work term (\f$\nabla p \cdot v\f$))
- * is
- \f{align*}{
- \frac{\partial}{\partial t} (\rho u) &= -\nabla \cdot (\rho v h) + \nabla \cdot (\lambda \nabla T)
- \f}
- *
- * Incompressible energy formulation is
- \f{align*}{
- \frac{\partial}{\partial t} (\rho u) = -\nabla \cdot (\rho v u) + \nabla \cdot (\lambda \nabla T)
- \f}
- *
- * \param flux The flux
- * \param fluxVars The flux variables.
- * \param phaseIdx The phase index
- */
- static void heatConvectionFlux(NumEqVector& flux,
- FluxVariables& fluxVars,
- int phaseIdx)
- {
- // internal energy used instead of enthalpy for incompressible flow
- auto upwindTerm = [phaseIdx](const auto& volVars)
- { return volVars.density(phaseIdx)*volVars.mobility(phaseIdx)*volVars.internalEnergy(phaseIdx); };
-
- flux[energyEqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
- }
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/nonisothermal/model.hh b/dumux/porousmediumflow/nonisothermal/model.hh
index 18eadc2149..828a1eccbc 100644
--- a/dumux/porousmediumflow/nonisothermal/model.hh
+++ b/dumux/porousmediumflow/nonisothermal/model.hh
@@ -19,13 +19,13 @@
* results in one energy conservation equation for the porous solid
* matrix and the fluids,
\f{align*}{
- \frac{\partial (\sum_\alpha \phi \varrho_\alpha u_\alpha S_\alpha )}{\partial t}
+ \phi \frac{\partial \sum_\alpha \varrho_\alpha S_\alpha (u_\alpha - \mathbf{g}\cdot\mathbf{x})}{\partial t}
& +
\frac{\partial \left((\left( 1 - \phi \right) \varrho_s c_s T\right)}{\partial t}
-
\sum_\alpha \nabla \cdot
\left\{
- \varrho_\alpha h_\alpha
+ \varrho_\alpha (h_\alpha - \mathbf{g}\cdot\mathbf{x})
\frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K}
\left( \nabla p_\alpha - \varrho_\alpha \mathbf{g} \right)
\right\} \\
@@ -49,6 +49,11 @@
* * \f$ \mathbf{g} \f$ is the gravitational acceleration vector,
* * \f$ q^h \f$ is a source or sink term.
*
+ * The gravitational potential is approximated by \f$\psi \approx \mathbf{g}\cdot\mathbf{x}\f$
+ * with \f$ \mathbf{g} = -\nabla\psi\f$ and where \f$\mathbf{x}\f$ is the position in space.
+ * The approximation is exact in case the gravitation vector is constant in space.
+ * The implementation via the gravitational potential in the formulation above also assumes
+ * that the gravitional potential is independent of time. This formulation is incorrect otherwise.
*/
#ifndef DUMUX_NONISOTHERMAL_MODEL_HH
diff --git a/dumux/porousmediumflow/richards/localresidual.hh b/dumux/porousmediumflow/richards/localresidual.hh
index f52fd97e1a..dbfd6d5346 100644
--- a/dumux/porousmediumflow/richards/localresidual.hh
+++ b/dumux/porousmediumflow/richards/localresidual.hh
@@ -108,8 +108,8 @@ public:
* volVars.saturation(liquidPhaseIdx);
//! The energy storage in the water, air and solid phase
- EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, liquidPhaseIdx);
- EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, gasPhaseIdx);
+ EnergyLocalResidual::fluidPhaseStorage(storage, problem, scv, volVars, liquidPhaseIdx);
+ EnergyLocalResidual::fluidPhaseStorage(storage, problem, scv, volVars, gasPhaseIdx);
EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
return storage;
--
GitLab