Commit 96ebfc73 authored by Katharina Heck's avatar Katharina Heck Committed by Kilian Weishaupt
Browse files

[porousmedium][compositional] add check for reference system in

localresidual and change units of the diffusive flux accoridingly. For
mass averaged systems the unit of the diffusive flux is in kg/s so has
to be devided by the molar mass for useMoles = true
parent 99d7595a
......@@ -137,6 +137,7 @@ public:
{
FluxVariables fluxVars;
fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
auto referenceSystemFormulation = FluxVariables::MolecularDiffusionType::referenceSystemFormulation();
// get upwind weights into local scope
NumEqVector flux(0.0);
......@@ -164,15 +165,44 @@ public:
const auto diffusionFluxesWPhase = fluxVars.molecularDiffusionFlux(wPhaseIdx);
// diffusive fluxes
Scalar jGW = diffusionFluxesWPhase[gCompIdx];
Scalar jNW = diffusionFluxesWPhase[nCompIdx];
Scalar jWW = -(jGW+jNW);
Scalar jGW = 0.0;
Scalar jNW = 0.0;
Scalar jWW = 0.0;
//check for the reference system and adapt units of the diffusive flux accordingly.
if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
{
jGW += diffusionFluxesWPhase[gCompIdx]/FluidSystem::molarMass(gCompIdx);
jNW += diffusionFluxesWPhase[nCompIdx]/FluidSystem::molarMass(nCompIdx);
}
else
{
jGW += diffusionFluxesWPhase[gCompIdx];
jNW += diffusionFluxesWPhase[nCompIdx];
}
jWW += -(jGW+jNW);
const auto diffusionFluxesGPhase = fluxVars.molecularDiffusionFlux(gPhaseIdx);
Scalar jWG = diffusionFluxesGPhase[wCompIdx];
Scalar jNG = diffusionFluxesGPhase[nCompIdx];
Scalar jGG = -(jWG+jNG);
Scalar jWG = 0.0;
Scalar jNG = 0.0;
Scalar jGG = 0.0;
//check for the reference system and adapt units of the diffusive flux accordingly.
if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
{
jWG += diffusionFluxesGPhase[wCompIdx]/FluidSystem::molarMass(wCompIdx);
jNG += diffusionFluxesGPhase[nCompIdx]/FluidSystem::molarMass(nCompIdx);
}
else
{
jWG += diffusionFluxesGPhase[wCompIdx];
jNG += diffusionFluxesGPhase[nCompIdx];
}
jGG += -(jWG+jNG);
// At the moment we do not consider diffusion in the NAPL phase
Scalar jWN = 0.0;
......
......@@ -174,19 +174,40 @@ public:
const auto diffusionFluxesWPhase = fluxVars.molecularDiffusionFlux(wPhaseIdx);
Scalar jNW = 0.0;
Scalar jWW = 0.0;
// diffusive fluxes
Scalar jNW = diffusionFluxesWPhase[nCompIdx];
Scalar jWW = -(jNW);
//check for the reference system and adapt units of the diffusive flux accordingly.
if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
jNW += diffusionFluxesWPhase[nCompIdx]/FluidSystem::molarMass(nCompIdx);
else
jNW += diffusionFluxesWPhase[nCompIdx];
jWW += -(jNW);
const auto diffusionFluxesGPhase = fluxVars.molecularDiffusionFlux(gPhaseIdx);
Scalar jWG = diffusionFluxesGPhase[wCompIdx];
Scalar jNG = -(jWG);
Scalar jWG = 0.0;
Scalar jNG = 0.0;
if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
jWG += diffusionFluxesGPhase[wCompIdx]/FluidSystem::molarMass(wCompIdx);
else
jWG += diffusionFluxesGPhase[wCompIdx];
jNG += -(jWG);
const auto diffusionFluxesNPhase = fluxVars.molecularDiffusionFlux(nPhaseIdx);
// At the moment we do not consider diffusion in the NAPL phase
Scalar jWN = diffusionFluxesNPhase[wCompIdx];
Scalar jNN = -jWN;
Scalar jWN = 0.0;
Scalar jNN = 0.0;
if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
jWN += diffusionFluxesNPhase[wCompIdx]/FluidSystem::molarMass(wCompIdx);
else
jWN += diffusionFluxesNPhase[wCompIdx];
jNN += -jWN;
flux[conti0EqIdx] += jWW+jWG+jWN;
flux[conti1EqIdx] += jNW+jNG+jNN;
......
......@@ -143,6 +143,7 @@ public:
{
FluxVariables fluxVars;
fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
auto referenceSystemFormulation = FluxVariables::MolecularDiffusionType::referenceSystemFormulation();
// get upwind weights into local scope
NumEqVector flux(0.0);
......@@ -170,8 +171,15 @@ public:
// diffusive fluxes (only for the component balances)
if(eqIdx != replaceCompEqIdx)
flux[eqIdx] += useMoles ? diffusiveFluxes[compIdx]
{
//check for the reference system and adapt units of the diffusive flux accordingly.
if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
flux[eqIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx)
: diffusiveFluxes[compIdx];
else
flux[eqIdx] += useMoles ? diffusiveFluxes[compIdx]
: diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
}
}
// in case one balance is substituted by the total mole balance
......@@ -184,8 +192,14 @@ public:
flux[replaceCompEqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
for(int compIdx = 0; compIdx < numComponents; ++compIdx)
flux[replaceCompEqIdx] += useMoles ? diffusiveFluxes[compIdx]
: diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
{
//check for the reference system and adapt units of the diffusive flux accordingly.
if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
flux[replaceCompEqIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx) : diffusiveFluxes[compIdx];
else
flux[replaceCompEqIdx] += useMoles ? diffusiveFluxes[compIdx]
: diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
}
}
//! Add advective phase energy fluxes. For isothermal model the contribution is zero.
......
......@@ -59,6 +59,7 @@ class NonEquilibriumLocalResidualImplementation<TypeTag, false>: public GetPropT
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using EnergyLocalResidual = GetPropType<TypeTag, Properties::EnergyLocalResidual>;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Indices = typename ModelTraits::Indices;
static constexpr int numPhases = ModelTraits::numFluidPhases();
......@@ -116,6 +117,10 @@ public:
flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
// diffusive fluxes (only for the component balances)
//check for the reference system and adapt units of the diffusive flux accordingly.
if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
flux[eqIdx] += diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx);
else
flux[eqIdx] += diffusiveFluxes[compIdx];
}
......@@ -271,7 +276,11 @@ public:
// do not add diffusive flux of main component, as that is not done in master as well
if (compIdx == phaseIdx)
continue;
flux[eqIdx] += diffusiveFluxes[compIdx];
//check for the reference system and adapt units of the diffusive flux accordingly.
if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
flux[eqIdx] += diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx);
else
flux[eqIdx] += diffusiveFluxes[compIdx];
}
//! Add advective phase energy fluxes. For isothermal model the contribution is zero.
EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
......
......@@ -128,7 +128,12 @@ public:
enthalpy += insideVolVars.enthalpy(phaseIdx);
else
enthalpy += outsideVolVars.enthalpy(phaseIdx);
flux[energyEq0Idx] += diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx)*enthalpy;
//check for the reference system and adapt units of the diffusive flux accordingly.
if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
flux[energyEq0Idx] += diffusiveFluxes[compIdx]*enthalpy;
else
flux[energyEq0Idx] += diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx)*enthalpy;
}
}
......
......@@ -148,7 +148,13 @@ public:
// for extended Richards we consider water vapor diffusion in air
if (enableWaterDiffusionInAir)
flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx]*FluidSystem::molarMass(liquidCompIdx);
{
//check for the reference system and adapt units of the diffusive flux accordingly.
if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx];
else
flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx]*FluidSystem::molarMass(liquidCompIdx);
}
//! Add advective phase energy fluxes for the water phase only. For isothermal model the contribution is zero.
EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, liquidPhaseIdx);
......
......@@ -138,7 +138,10 @@ public:
// advective fluxes
flux[compIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
// diffusive fluxes
flux[compIdx] += diffusiveFluxes[compIdx];
if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
flux[compIdx] += diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx);
else
flux[compIdx] += diffusiveFluxes[compIdx];
}
}
// formulation with mass balances
......@@ -153,7 +156,10 @@ public:
// advective fluxes
flux[compIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
// diffusive fluxes
flux[compIdx] += diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
flux[compIdx] += diffusiveFluxes[compIdx];
else
flux[compIdx] += diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
}
}
......@@ -237,14 +243,28 @@ public:
const auto advDerivIJ = volFlux*rho(outsideVolVars)*outsideWeight;
// diffusive term
auto referenceSystemFormulation = FluxVariables::MolecularDiffusionType::referenceSystemFormulation();
const auto& fluxCache = elemFluxVarsCache[scvf];
const auto rhoMolar = 0.5*(insideVolVars.molarDensity() + outsideVolVars.molarDensity());
const auto rhoInside = (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged) ? insideVolVars.density(phaseIdx) : insideVolVars.molarDensity(phaseIdx);
const auto rhoOutside = (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged) ? outsideVolVars.density(phaseIdx) : outsideVolVars.molarDensity(phaseIdx);
const auto massOrMolarDensity = 0.5*(rhoInside + rhoOutside);
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
// diffusive term
const auto diffDeriv = useMoles ? rhoMolar*fluxCache.diffusionTij(phaseIdx, compIdx)
: rhoMolar*fluxCache.diffusionTij(phaseIdx, compIdx)*FluidSystem::molarMass(compIdx);
auto diffDeriv = 0.0;
if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
{
diffDeriv = useMoles ? massOrMolarDensity*fluxCache.diffusionTij(phaseIdx, compIdx)/FluidSystem::molarMass(compIdx)
: massOrMolarDensity*fluxCache.diffusionTij(phaseIdx, compIdx);
}
else
{
diffDeriv = useMoles ? massOrMolarDensity*fluxCache.diffusionTij(phaseIdx, compIdx)
: massOrMolarDensity*fluxCache.diffusionTij(phaseIdx, compIdx)*FluidSystem::molarMass(compIdx);
}
derivativeMatrices[scvf.insideScvIdx()][compIdx][compIdx] += (advDerivII + diffDeriv);
derivativeMatrices[scvf.outsideScvIdx()][compIdx][compIdx] += (advDerivIJ - diffDeriv);
......@@ -282,6 +302,7 @@ public:
const auto advDerivIJ = volFlux*rho(outsideVolVars)*outsideWeight;
// diffusive term
auto referenceSystemFormulation = FluxVariables::MolecularDiffusionType::referenceSystemFormulation();
using DiffusionType = GetPropType<T, Properties::MolecularDiffusionType>;
const auto ti = DiffusionType::calculateTransmissibilities(problem,
element,
......@@ -298,8 +319,13 @@ public:
for (const auto& scv : scvs(fvGeometry))
{
// diffusive term
const auto diffDeriv = useMoles ? ti[compIdx][scv.indexInElement()]
: ti[compIdx][scv.indexInElement()]*FluidSystem::molarMass(compIdx);
auto diffDeriv = 0.0;
if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
diffDeriv += useMoles ? ti[compIdx][scv.indexInElement()]/FluidSystem::molarMass(compIdx)
: ti[compIdx][scv.indexInElement()];
else
diffDeriv += useMoles ? ti[compIdx][scv.indexInElement()]
: ti[compIdx][scv.indexInElement()]*FluidSystem::molarMass(compIdx);
A[insideScv.dofIndex()][scv.dofIndex()][compIdx][compIdx] += diffDeriv;
A[outsideScv.dofIndex()][scv.dofIndex()][compIdx][compIdx] -= diffDeriv;
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment