From 70bb0704fe274d6cd4b5f9260e8dc80f10f1f689 Mon Sep 17 00:00:00 2001
From: Klaus Mosthaf <klmos@env.dtu.dk>
Date: Mon, 11 Jul 2011 09:27:31 +0000
Subject: [PATCH] Added the possibility to replace one of the two component
mass balances by a total mass balance. This was done by adding a property
ReplaceCompIdx, which defines which equation is replaced. NOTE: adjust the
boundary conditions and source/sink values if one component balance is
replaced. If ReplaceCompIdx is >= numComponents (set as default), the model
does the same as before. This change was neccessary for the model coupling of
2cnistokes (uses a total mass balance) and 2p2cni.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6195 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/2p2c/2p2clocalresidual.hh | 82 +++++++++++++++++---
dumux/boxmodels/2p2c/2p2cproperties.hh | 1 +
dumux/boxmodels/2p2c/2p2cpropertydefaults.hh | 3 +
3 files changed, 77 insertions(+), 9 deletions(-)
diff --git a/dumux/boxmodels/2p2c/2p2clocalresidual.hh b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
index f7c4486031..136b3a2642 100644
--- a/dumux/boxmodels/2p2c/2p2clocalresidual.hh
+++ b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
@@ -33,11 +33,8 @@
#include <dumux/common/math.hh>
#include "2p2cproperties.hh"
-
#include "2p2cvolumevariables.hh"
-
#include "2p2cfluxvariables.hh"
-
#include "2p2cnewtoncontroller.hh"
#include <iostream>
@@ -121,6 +118,8 @@ protected:
static constexpr Scalar mobilityUpwindAlpha =
GET_PROP_VALUE(TypeTag, PTAG(MobilityUpwindAlpha));
+ static constexpr unsigned int replaceCompEqIdx =
+ GET_PROP_VALUE(TypeTag, PTAG(ReplaceCompEqIdx));
public:
/*!
@@ -174,15 +173,21 @@ public:
// compute storage term of all components within all phases
result = 0;
+
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ for (int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx)
{
int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
result[eqIdx] += volVars.density(phaseIdx)
* volVars.saturation(phaseIdx)
* volVars.fluidState().massFrac(phaseIdx, compIdx);
}
+ // this is only processed, if one component mass balance equation
+ // is replaced by a total mass balance equation
+ if (replaceCompEqIdx < numComponents)
+ result[replaceCompEqIdx] += volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx);
}
result *= volVars.porosity();
}
@@ -230,7 +235,7 @@ public:
const VolumeVariables &dn =
this->curVolVars_(vars.downstreamIdx(phaseIdx));
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ for (int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx)
{
int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
// add advective flux of current component in current
@@ -247,6 +252,7 @@ public:
- mobilityUpwindAlpha) * (dn.density(phaseIdx)
* dn.mobility(phaseIdx) * dn.fluidState().massFrac(
phaseIdx, compIdx));
+
Valgrind::CheckDefined(vars.KmvpNormal(phaseIdx));
Valgrind::CheckDefined(up.density(phaseIdx));
Valgrind::CheckDefined(up.mobility(phaseIdx));
@@ -255,6 +261,29 @@ public:
Valgrind::CheckDefined(dn.mobility(phaseIdx));
Valgrind::CheckDefined(dn.fluidState().massFrac(phaseIdx, compIdx));
}
+ // flux of the total mass balance;
+ // this is only processed, if one component mass balance equation
+ // is replaced by a total mass balance equation
+ if (replaceCompEqIdx < numComponents)
+ {
+ // upstream vertex
+ if (mobilityUpwindAlpha > 0.0)
+ flux[replaceCompEqIdx] += vars.KmvpNormal(phaseIdx)
+ * mobilityUpwindAlpha *
+ (up.density(phaseIdx) * up.mobility(phaseIdx));
+ // downstream vertex
+ if (mobilityUpwindAlpha < 1.0)
+ flux[replaceCompEqIdx] += vars.KmvpNormal(phaseIdx) * (1
+ - mobilityUpwindAlpha) * (dn.density(phaseIdx)
+ * dn.mobility(phaseIdx));
+ Valgrind::CheckDefined(vars.KmvpNormal(phaseIdx));
+ Valgrind::CheckDefined(up.density(phaseIdx));
+ Valgrind::CheckDefined(up.mobility(phaseIdx));
+ Valgrind::CheckDefined(dn.density(phaseIdx));
+ Valgrind::CheckDefined(dn.mobility(phaseIdx));
+
+ }
+
}
}
@@ -272,16 +301,22 @@ public:
- vars.porousDiffCoeff(lPhaseIdx) *
vars.molarDensityAtIP(lPhaseIdx) *
(vars.molarConcGrad(lPhaseIdx) * vars.face().normal);
- flux[contiGEqIdx] += tmp * FluidSystem::molarMass(gCompIdx);
- flux[contiLEqIdx] -= tmp * FluidSystem::molarMass(lCompIdx);
+ // add the diffusive fluxes only to the component mass balance
+ if (replaceCompEqIdx != contiGEqIdx)
+ flux[contiGEqIdx] += tmp * FluidSystem::molarMass(gCompIdx);
+ if (replaceCompEqIdx != contiLEqIdx)
+ flux[contiLEqIdx] -= tmp * FluidSystem::molarMass(lCompIdx);
// add diffusive flux of liquid component in gas phase
tmp =
- vars.porousDiffCoeff(gPhaseIdx) *
vars.molarDensityAtIP(gPhaseIdx) *
(vars.molarConcGrad(gPhaseIdx) * vars.face().normal);
- flux[contiLEqIdx] += tmp * FluidSystem::molarMass(lCompIdx);
- flux[contiGEqIdx] -= tmp * FluidSystem::molarMass(gCompIdx);
+ // add the diffusive fluxes only to the component mass balance
+ if (replaceCompEqIdx != contiLEqIdx)
+ flux[contiLEqIdx] += tmp * FluidSystem::molarMass(lCompIdx);
+ if (replaceCompEqIdx != contiGEqIdx)
+ flux[contiGEqIdx] -= tmp * FluidSystem::molarMass(gCompIdx);
}
/*!
@@ -323,6 +358,35 @@ protected:
}
}
+ /*!
+ * \brief Return the equation index of the first mass balance equation
+ * of the component (used for loops); if one component mass balance
+ * is replaced by the total mass balance, this is the index
+ * of the remaining component mass balance equation.
+ */
+ constexpr unsigned int contiCompIdx1_() const {
+ switch (replaceCompEqIdx)
+ {
+ case contiLEqIdx: return contiGEqIdx;
+ case contiGEqIdx: return contiLEqIdx;
+ default: return 0;
+ }
+ }
+
+ /*!
+ * \brief Return the equation index of the second mass balance
+ * of the component (used for loops);
+ * if one component mass balance is replaced by the total mass balance
+ * (replaceCompEqIdx < 2), this index is the same as contiCompIdx1().
+ */
+ constexpr unsigned int contiCompIdx2_() const {
+ switch (replaceCompEqIdx)
+ {
+ case contiLEqIdx: return contiGEqIdx;
+ case contiGEqIdx: return contiLEqIdx;
+ default: return numComponents-1;
+ }
+ }
Implementation *asImp_()
{
diff --git a/dumux/boxmodels/2p2c/2p2cproperties.hh b/dumux/boxmodels/2p2c/2p2cproperties.hh
index cd6436d5e2..18821bf6fa 100644
--- a/dumux/boxmodels/2p2c/2p2cproperties.hh
+++ b/dumux/boxmodels/2p2c/2p2cproperties.hh
@@ -61,6 +61,7 @@ NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extrac
NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
NEW_PROP_TAG(MobilityUpwindAlpha); //!< The value of the upwind parameter for the mobility
+NEW_PROP_TAG(ReplaceCompEqIdx); //!< The index of the total mass balance equation, if one component balance is replaced (ReplaceCompEqIdx < NumComponents)
}
}
diff --git a/dumux/boxmodels/2p2c/2p2cpropertydefaults.hh b/dumux/boxmodels/2p2c/2p2cpropertydefaults.hh
index 8f7aa9e8c6..623003a867 100644
--- a/dumux/boxmodels/2p2c/2p2cpropertydefaults.hh
+++ b/dumux/boxmodels/2p2c/2p2cpropertydefaults.hh
@@ -91,6 +91,9 @@ SET_INT_PROP(BoxTwoPTwoC,
Formulation,
TwoPTwoCFormulation::plSg);
+//! set as default that no component mass balance is replaced by the total mass balance
+SET_INT_PROP(BoxTwoPTwoC, ReplaceCompEqIdx, 2);
+
/*!
* \brief Set the property for the material law by retrieving it from
* the spatial parameters.
--
GitLab