From 3556eac76c87ba20c23d412ab99319bb3ff47aff Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Mon, 19 Nov 2018 09:16:51 +0100
Subject: [PATCH] [mpnc][test] use extra combustionlocalresidual for combustion
test as this is has some unique qboil functions. Cleanup normal
nonequilibrium thermal localresidual, does now only normal energy exchange
between phases
---
.../nonequilibrium/thermal/localresidual.hh | 120 +---------
.../thermalnonequilibrium/CMakeLists.txt | 1 +
.../combustionlocalresidual.hh | 225 ++++++++++++++++++
.../implicit/thermalnonequilibrium/problem.hh | 5 +
4 files changed, 238 insertions(+), 113 deletions(-)
create mode 100644 test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh
diff --git a/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh b/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
index f714851509..4047df849d 100644
--- a/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
+++ b/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
@@ -66,8 +66,6 @@ class EnergyLocalResidualNonEquilibrium<TypeTag, 1/*numEnergyEqFluid*/>
enum { numPhases = ModelTraits::numPhases() };
enum { numComponents = ModelTraits::numComponents() };
- enum { phase0Idx = FluidSystem::phase0Idx};
- enum { phase1Idx = FluidSystem::phase1Idx};
public:
//! The energy storage in the fluid phase with index phaseIdx
@@ -147,12 +145,12 @@ public:
FluxVariables& fluxVars)
{
//in case we have one energy equation for more than one fluid phase we use an effective law in the nonequilibrium fourierslaw
- flux[energyEq0Idx] += fluxVars.heatConductionFlux(0);
+ flux[energyEq0Idx] += fluxVars.heatConductionFlux(0);
//heat conduction for the solid phases
- for(int sPhaseIdx=0; sPhaseIdx<numEnergyEqSolid; ++sPhaseIdx)
- {
- flux[energyEqSolidIdx+sPhaseIdx] += fluxVars.heatConductionFlux(numPhases + sPhaseIdx);
- }
+ for(int sPhaseIdx=0; sPhaseIdx<numEnergyEqSolid; ++sPhaseIdx)
+ {
+ flux[energyEqSolidIdx+sPhaseIdx] += fluxVars.heatConductionFlux(numPhases + sPhaseIdx);
+ }
}
/*!
@@ -168,7 +166,6 @@ public:
{
//specialization for 2 fluid phases
const auto& volVars = elemVolVars[scv];
- const auto& fs = volVars.fluidState() ;
const Scalar characteristicLength = volVars.characteristicLength() ;
//interfacial area
@@ -179,71 +176,15 @@ public:
const Scalar TFluid = volVars.temperatureFluid(0);
const Scalar TSolid = volVars.temperatureSolid();
- const Scalar satW = fs.saturation(phase0Idx) ;
- const Scalar satN = fs.saturation(phase1Idx) ;
-
- const Scalar eps = 1e-6 ;
Scalar solidToFluidEnergyExchange ;
- Scalar fluidConductivity ;
- if (satW < 1.0 - eps)
- fluidConductivity = volVars.fluidThermalConductivity(phase1Idx) ;
- else if (satW >= 1.0 - eps)
- fluidConductivity = volVars.fluidThermalConductivity(phase0Idx) ;
- else
- DUNE_THROW(Dune::NotImplemented,
- "wrong range");
+ const Scalar fluidConductivity = volVars.fluidThermalConductivity(0) ;
const Scalar factorEnergyTransfer = volVars.factorEnergyTransfer() ;
solidToFluidEnergyExchange = factorEnergyTransfer * (TSolid - TFluid) / characteristicLength * as * fluidConductivity ;
- const Scalar epsRegul = 1e-3 ;
-
- if (satW < (0 + eps) )
- {
- solidToFluidEnergyExchange *= volVars.nusseltNumber(phase1Idx) ;
- }
- else if ( (satW >= 0 + eps) and (satW < 1.0-eps) )
- {
- solidToFluidEnergyExchange *= (volVars.nusseltNumber(phase1Idx) * satN );
- Scalar qBoil ;
- if (satW<=epsRegul)
- {// regularize
- typedef Dumux::Spline<Scalar> Spline;
- Spline sp(0.0, epsRegul, // x1, x2
- QBoilFunc(volVars, 0.0), QBoilFunc(volVars, epsRegul), // y1, y2
- 0.0,dQBoil_dSw(volVars, epsRegul)); // m1, m2
-
- qBoil = sp.eval(satW) ;
- }
-
- else if (satW>= (1.0-epsRegul) )
- {// regularize
- typedef Dumux::Spline<Scalar> Spline;
- Spline sp(1.0-epsRegul, 1.0, // x1, x2
- QBoilFunc(volVars, 1.0-epsRegul), 0.0, // y1, y2
- dQBoil_dSw(volVars, 1.0-epsRegul), 0.0 ); // m1, m2
-
- qBoil = sp.eval(satW) ;
- }
- else
- {
- qBoil = QBoilFunc(volVars, satW) ;
- }
-
- solidToFluidEnergyExchange += qBoil;
- }
- else if (satW >= 1.0-eps)
- {
- solidToFluidEnergyExchange *= volVars.nusseltNumber(phase0Idx) ;
- }
- else
- DUNE_THROW(Dune::NotImplemented,
- "wrong range");
- using std::isfinite;
- if (!isfinite(solidToFluidEnergyExchange))
- DUNE_THROW(NumericalProblem, "Calculated non-finite source, " << "TFluid="<< TFluid << " TSolid="<< TSolid);
+ solidToFluidEnergyExchange *= volVars.nusseltNumber(0) ;
for(int energyEqIdx =0; energyEqIdx<numEnergyEqFluid+numEnergyEqSolid; ++energyEqIdx)
{
@@ -261,53 +202,6 @@ public:
} // end switch
}// end energyEqIdx
}// end source
-
-
- /*! \brief Calculate the energy transfer during boiling, i.e. latent heat
- *
- * \param volVars The volume variables
- * \param satW The wetting phase saturation. Not taken from volVars, because we regularize.
- */
- static Scalar QBoilFunc(const VolumeVariables & volVars,
- const Scalar satW)
- {
- // using saturation as input (instead of from volVars)
- // in order to make regularization (evaluation at different points) easyer
- const auto& fs = volVars.fluidState() ;
- const Scalar g( 9.81 ) ;
- const Scalar gamma(0.0589) ;
- const Scalar TSolid = volVars.temperatureSolid();
- const Scalar characteristicLength = volVars.characteristicLength() ;
- using std::pow;
- const Scalar as = 6.0 * (1.0-volVars.porosity()) / characteristicLength ;
- const Scalar mul = fs.viscosity(phase0Idx) ;
- const Scalar deltahv = fs.enthalpy(phase1Idx) - fs.enthalpy(phase0Idx);
- const Scalar deltaRho = fs.density(phase0Idx) - fs.density(phase1Idx) ;
- const Scalar firstBracket = pow(g * deltaRho / gamma, 0.5);
- const Scalar cp = FluidSystem::heatCapacity(fs, phase0Idx) ;
- // This use of Tsat is only justified if the fluid is always boiling (tsat equals boiling conditions)
- // If a different state is to be simulated, please use the actual fluid temperature instead.
- const Scalar Tsat = FluidSystem::vaporTemperature(fs, phase1Idx ) ;
- const Scalar deltaT = TSolid - Tsat ;
- const Scalar secondBracket = pow( (cp *deltaT / (0.006 * deltahv) ) , 3.0 ) ;
- const Scalar Prl = volVars.prandtlNumber(phase0Idx) ;
- const Scalar thirdBracket = pow( 1/Prl , (1.7/0.33) );
- const Scalar QBoil = satW * as * mul * deltahv * firstBracket * secondBracket * thirdBracket ;
- return QBoil;
- }
-
- /*! \brief Calculate the derivative of the energy transfer function during boiling. Needed for regularization.
- *
- * \param volVars The volume variables
- * \param satW The wetting phase saturation. Not taken from volVars, because we regularize.
- */
- static Scalar dQBoil_dSw(const VolumeVariables & volVars,
- const Scalar satW)
- {
- // on the fly derive w.r.t. Sw.
- // Only linearly depending on it (directly)
- return (QBoilFunc(volVars, satW) / satW ) ;
- }
};
template<class TypeTag>
diff --git a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/CMakeLists.txt b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/CMakeLists.txt
index 9fff88d75d..65f3ad9ac3 100644
--- a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/CMakeLists.txt
+++ b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/CMakeLists.txt
@@ -15,6 +15,7 @@ dune_add_test(COMPILE_ONLY # since it currently fails miserably with very differ
#install sources
install(FILES
combustionfluidsystem.hh
+combustionlocalresidual.hh
problem.hh
spatialparams.hh
main.cc
diff --git a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh
new file mode 100644
index 0000000000..4e2d8ff702
--- /dev/null
+++ b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh
@@ -0,0 +1,225 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup PorousmediumThermalNonEquilibriumModel
+ * \brief This file contains the parts of the local residual to
+ * calculate the heat conservation in the thermal non-equilibrium model.
+ */
+#ifndef DUMUX_ENERGY_COMBUSTION_LOCAL_RESIDUAL_HH
+#define DUMUX_ENERGY_COMBUSTION_LOCAL_RESIDUAL_HH
+
+#include <cmath>
+#include <dumux/common/spline.hh>
+#include <dumux/common/exceptions.hh>
+#include <dumux/common/properties.hh>
+
+#include <dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup PorousmediumThermalNonEquilibriumModel
+ * \brief This file contains the parts of the local residual to
+ * calculate the heat conservation in the thermal non-equilibrium model.
+ */
+template<class TypeTag>
+class CombustionEnergyLocalResidual
+: public EnergyLocalResidualNonEquilibrium<TypeTag, 1/*numEnergyEqFluid*/>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ using Indices = typename ModelTraits::Indices;
+
+ enum { numEnergyEqFluid = ModelTraits::numEnergyEqFluid() };
+ enum { numEnergyEqSolid = ModelTraits::numEnergyEqSolid() };
+ enum { energyEq0Idx = Indices::energyEq0Idx };
+ enum { energyEqSolidIdx = Indices::energyEqSolidIdx};
+
+ enum { numComponents = ModelTraits::numComponents() };
+
+public:
+ /*!
+ * \brief Calculate the source term of the equation
+ *
+ * \param scv The sub-control volume over which we integrate the source term
+ */
+ static void computeSourceEnergy(NumEqVector& source,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv)
+ {
+ //specialization for 2 fluid phases
+ const auto& volVars = elemVolVars[scv];
+ const auto& fs = volVars.fluidState() ;
+ const Scalar characteristicLength = volVars.characteristicLength() ;
+
+ //interfacial area
+ // Shi & Wang, Transport in porous media (2011)
+ const Scalar as = 6.0 * (1.0-volVars.porosity()) / characteristicLength ;
+
+ //temperature fluid is the same for both fluids
+ const Scalar TFluid = volVars.temperatureFluid(0);
+ const Scalar TSolid = volVars.temperatureSolid();
+
+ const Scalar satW = fs.saturation(0) ;
+ const Scalar satN = fs.saturation(1) ;
+
+ const Scalar eps = 1e-6 ;
+ Scalar solidToFluidEnergyExchange ;
+
+ Scalar fluidConductivity ;
+ if (satW < 1.0 - eps)
+ fluidConductivity = volVars.fluidThermalConductivity(1) ;
+ else if (satW >= 1.0 - eps)
+ fluidConductivity = volVars.fluidThermalConductivity(0) ;
+ else
+ DUNE_THROW(Dune::NotImplemented,
+ "wrong range");
+
+ const Scalar factorEnergyTransfer = volVars.factorEnergyTransfer() ;
+
+ solidToFluidEnergyExchange = factorEnergyTransfer * (TSolid - TFluid) / characteristicLength * as * fluidConductivity ;
+ const Scalar epsRegul = 1e-3 ;
+
+ if (satW < (0 + eps) )
+ {
+ solidToFluidEnergyExchange *= volVars.nusseltNumber(1) ;
+ }
+ else if ( (satW >= 0 + eps) and (satW < 1.0-eps) )
+ {
+ solidToFluidEnergyExchange *= (volVars.nusseltNumber(1) * satN );
+ Scalar qBoil ;
+ if (satW<=epsRegul)
+ {// regularize
+ typedef Dumux::Spline<Scalar> Spline;
+ Spline sp(0.0, epsRegul, // x1, x2
+ QBoilFunc(volVars, 0.0), QBoilFunc(volVars, epsRegul), // y1, y2
+ 0.0,dQBoil_dSw(volVars, epsRegul)); // m1, m2
+
+ qBoil = sp.eval(satW) ;
+ }
+
+ else if (satW>= (1.0-epsRegul) )
+ {// regularize
+ typedef Dumux::Spline<Scalar> Spline;
+ Spline sp(1.0-epsRegul, 1.0, // x1, x2
+ QBoilFunc(volVars, 1.0-epsRegul), 0.0, // y1, y2
+ dQBoil_dSw(volVars, 1.0-epsRegul), 0.0 ); // m1, m2
+
+ qBoil = sp.eval(satW) ;
+ }
+ else
+ {
+ qBoil = QBoilFunc(volVars, satW) ;
+ }
+
+ solidToFluidEnergyExchange += qBoil;
+ }
+ else if (satW >= 1.0-eps)
+ {
+ solidToFluidEnergyExchange *= volVars.nusseltNumber(0) ;
+ }
+ else
+ DUNE_THROW(Dune::NotImplemented,
+ "wrong range");
+
+ using std::isfinite;
+ if (!isfinite(solidToFluidEnergyExchange))
+ DUNE_THROW(NumericalProblem, "Calculated non-finite source, " << "TFluid="<< TFluid << " TSolid="<< TSolid);
+
+ for(int energyEqIdx =0; energyEqIdx<numEnergyEqFluid+numEnergyEqSolid; ++energyEqIdx)
+ {
+ switch (energyEqIdx)
+ {
+ case 0 :
+ source[energyEq0Idx + energyEqIdx] += solidToFluidEnergyExchange;
+ break;
+ case 1 :
+ source[energyEq0Idx + energyEqIdx] -= solidToFluidEnergyExchange;
+ break;
+ default:
+ DUNE_THROW(Dune::NotImplemented,
+ "wrong index");
+ } // end switch
+ }// end energyEqIdx
+ }// end source
+
+
+ /*! \brief Calculate the energy transfer during boiling, i.e. latent heat
+ *
+ * \param volVars The volume variables
+ * \param satW The wetting phase saturation. Not taken from volVars, because we regularize.
+ */
+ static Scalar QBoilFunc(const VolumeVariables & volVars,
+ const Scalar satW)
+ {
+ // using saturation as input (instead of from volVars)
+ // in order to make regularization (evaluation at different points) easyer
+ const auto& fs = volVars.fluidState() ;
+ const Scalar g( 9.81 ) ;
+ const Scalar gamma(0.0589) ;
+ const Scalar TSolid = volVars.temperatureSolid();
+ const Scalar characteristicLength = volVars.characteristicLength() ;
+ using std::pow;
+ const Scalar as = 6.0 * (1.0-volVars.porosity()) / characteristicLength ;
+ const Scalar mul = fs.viscosity(0) ;
+ const Scalar deltahv = fs.enthalpy(1) - fs.enthalpy(0);
+ const Scalar deltaRho = fs.density(0) - fs.density(1) ;
+ const Scalar firstBracket = pow(g * deltaRho / gamma, 0.5);
+ const Scalar cp = FluidSystem::heatCapacity(fs, 0) ;
+ // This use of Tsat is only justified if the fluid is always boiling (tsat equals boiling conditions)
+ // If a different state is to be simulated, please use the actual fluid temperature instead.
+ const Scalar Tsat = FluidSystem::vaporTemperature(fs, 1 ) ;
+ const Scalar deltaT = TSolid - Tsat ;
+ const Scalar secondBracket = pow( (cp *deltaT / (0.006 * deltahv) ) , 3.0 ) ;
+ const Scalar Prl = volVars.prandtlNumber(0) ;
+ const Scalar thirdBracket = pow( 1/Prl , (1.7/0.33) );
+ const Scalar QBoil = satW * as * mul * deltahv * firstBracket * secondBracket * thirdBracket ;
+ return QBoil;
+ }
+
+ /*! \brief Calculate the derivative of the energy transfer function during boiling. Needed for regularization.
+ *
+ * \param volVars The volume variables
+ * \param satW The wetting phase saturation. Not taken from volVars, because we regularize.
+ */
+ static Scalar dQBoil_dSw(const VolumeVariables & volVars,
+ const Scalar satW)
+ {
+ // on the fly derive w.r.t. Sw.
+ // Only linearly depending on it (directly)
+ return (QBoilFunc(volVars, satW) / satW ) ;
+ }
+};
+} // end namespace Dumux
+
+#endif //
diff --git a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh
index 596bd1b626..2ef8a197bd 100644
--- a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh
+++ b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh
@@ -44,6 +44,7 @@
#include "spatialparams.hh"
#include "combustionfluidsystem.hh"
+#include "combustionlocalresidual.hh"
namespace Dumux {
@@ -154,6 +155,10 @@ private:
public:
using type = CompositionalSolidState<Scalar, SolidSystem>;
};
+
+template<class TypeTag>
+struct EnergyLocalResidual<TypeTag, TTag::CombustionOneComponent>
+{ using type = CombustionEnergyLocalResidual<TypeTag>; };
}
/*!
* \ingroup MPNCTests
--
GitLab