diff --git a/dumux/implicit/mpnc/diffusion/diffusion.hh b/dumux/implicit/mpnc/diffusion/diffusion.hh index f61e8108d55ff5ad802ba149df63a1b74bc9f7c5..5d470a10be04df0c121132fd3b9baf3311d7f9fd 100644 --- a/dumux/implicit/mpnc/diffusion/diffusion.hh +++ b/dumux/implicit/mpnc/diffusion/diffusion.hh @@ -43,6 +43,8 @@ class MPNCDiffusion enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents)}; enum { nPhaseIdx = FluidSystem::nPhaseIdx }; enum { wPhaseIdx = FluidSystem::wPhaseIdx }; + enum { wCompIdx = FluidSystem::wCompIdx }; + enum { nCompIdx = FluidSystem::nCompIdx }; typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables; typedef Dune::FieldMatrix<Scalar, numComponents, numComponents> ComponentMatrix; @@ -52,7 +54,7 @@ public: static void flux(ComponentVector &fluxes, const unsigned int phaseIdx, const FluxVariables &fluxVars, - const Scalar molarDensity) + const Scalar molarDensity ) { if (phaseIdx == nPhaseIdx) gasFlux_(fluxes, fluxVars, molarDensity); @@ -86,35 +88,48 @@ protected: const FluxVariables &fluxVars, const Scalar molarDensity) { - // Stefan-Maxwell equation - // - // See: R. Reid, et al.: "The Properties of Liquids and - // Gases", 4th edition, 1987, McGraw-Hill, p 596 - - // TODO: tensorial diffusion coefficients - ComponentMatrix M(0); - - for (int compIIdx = 0; compIIdx < numComponents - 1; ++compIIdx) { - for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) { - Scalar Dij = fluxVars.porousDiffCoeffG(compIIdx, compJIdx); - if (Dij) { - M[compIIdx][compJIdx] += fluxVars.moleFraction(nPhaseIdx, compIIdx) / Dij; - M[compIIdx][compIIdx] -= fluxVars.moleFraction(nPhaseIdx, compJIdx) / Dij; + // Alternative: use Fick Diffusion: no mutual influence of diffusion + if (GET_PROP_VALUE(TypeTag, UseMaxwellDiffusion) ){ + // Stefan-Maxwell equation + // + // See: R. Reid, et al.: "The Properties of Liquids and + // Gases", 4th edition, 1987, McGraw-Hill, p 596 + + // TODO: tensorial diffusion coefficients + ComponentMatrix M(0); + + for (int compIIdx = 0; compIIdx < numComponents - 1; ++compIIdx) { + for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) { + Scalar Dij = fluxVars.porousDiffCoeffG(compIIdx, compJIdx); + if (Dij) { + M[compIIdx][compJIdx] += fluxVars.moleFraction(nPhaseIdx, compIIdx) / Dij; + M[compIIdx][compIIdx] -= fluxVars.moleFraction(nPhaseIdx, compJIdx) / Dij; + } } } - } - for (int compIIdx = 0; compIIdx < numComponents; ++compIIdx) { - M[numComponents - 1][compIIdx] = 1.0; - } + for (int compIIdx = 0; compIIdx < numComponents; ++compIIdx) { + M[numComponents - 1][compIIdx] = 1.0; + } - ComponentVector rightHandSide ; // see source cited above - for (int compIIdx = 0; compIIdx < numComponents - 1; ++compIIdx) { - rightHandSide[compIIdx] = molarDensity*(fluxVars.moleFractionGrad(nPhaseIdx, compIIdx)*fluxVars.face().normal); - } - rightHandSide[numComponents - 1] = 0.0; + ComponentVector rightHandSide ; // see source cited above + for (int compIIdx = 0; compIIdx < numComponents - 1; ++compIIdx) { + rightHandSide[compIIdx] = molarDensity*(fluxVars.moleFractionGrad(nPhaseIdx, compIIdx)*fluxVars.face().normal); + } + rightHandSide[numComponents - 1] = 0.0; - M.solve(fluxes, rightHandSide); + M.solve(fluxes, rightHandSide); + } + else{// Fick Diffusion + for (int compIdx = 0; compIdx < numComponents; ++compIdx) { + // TODO: tensorial diffusion coefficients + const Scalar xGrad = fluxVars.moleFractionGrad(nPhaseIdx, compIdx)*fluxVars.face().normal; + fluxes[compIdx] = + - xGrad * + molarDensity + * fluxVars.porousDiffCoeffG(compIdx, nCompIdx) ; // this is == 0 for nComp==comp, i.e. no diffusion of the main component of the phase + } + } } // return whether a concentration can be assumed to be a trace diff --git a/dumux/implicit/mpnc/mpncproperties.hh b/dumux/implicit/mpnc/mpncproperties.hh index 5292d1c40831a91266bb91b352b104e7919beb83..de3c264dc0b054a21f34734eb6e1a3e0e39b9095 100644 --- a/dumux/implicit/mpnc/mpncproperties.hh +++ b/dumux/implicit/mpnc/mpncproperties.hh @@ -113,6 +113,9 @@ NEW_PROP_TAG(EnableKinetic); //! Enable kinetic resolution of energy transfer processes? NEW_PROP_TAG(EnableKineticEnergy); +//! Enable Maxwell Diffusion? (If false: use Fickian Diffusion) Maxwell incorporated the mutual influences of multiple diffusing components. However, Fick seems to be more robust. +NEW_PROP_TAG(UseMaxwellDiffusion); + //! Enable gravity? NEW_PROP_TAG(ProblemEnableGravity); diff --git a/dumux/implicit/mpnc/mpncpropertydefaults.hh b/dumux/implicit/mpnc/mpncpropertydefaults.hh index 96d4eb6b6fa8cec78e290d2319ef2d5ebc1023c2..9a7e1700a4fa923399a057fce6e4cc736db4415c 100644 --- a/dumux/implicit/mpnc/mpncpropertydefaults.hh +++ b/dumux/implicit/mpnc/mpncpropertydefaults.hh @@ -144,6 +144,9 @@ SET_BOOL_PROP(MPNC, EnableKinetic, false); //! disable kinetic energy transfer by default SET_BOOL_PROP(MPNC, EnableKineticEnergy, false); +//! disable Maxwell Diffusion by default: use Fick +SET_BOOL_PROP(MPNC, UseMaxwellDiffusion, false); + //! enable smooth upwinding by default SET_BOOL_PROP(MPNC, ImplicitEnableSmoothUpwinding, false);