diff --git a/appl/lecture/msm/1p2cvs2p/watercontaminantfluidsystem.hh b/appl/lecture/msm/1p2cvs2p/watercontaminantfluidsystem.hh index 91169019d66a3c955a778f778917f268285f535a..99f490b5ed9077f9a11257d2348802a955bc5d2f 100644 --- a/appl/lecture/msm/1p2cvs2p/watercontaminantfluidsystem.hh +++ b/appl/lecture/msm/1p2cvs2p/watercontaminantfluidsystem.hh @@ -282,7 +282,7 @@ public: /*! * \brief Given a phase's composition, temperature, pressure and - * density, calculate its specific internal energy [J/kg]. + * density, calculate its specific enthalpy [J/kg]. * * \todo This fluid system neglects the contribution of * gas-molecules in the liquid phase. This contribution is @@ -294,12 +294,12 @@ public: * \param phaseIdx The index of the fluid phase to consider */ template <class FluidState> - static Scalar internalEnergy(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + static Scalar enthalpy(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { // TODO! - DUNE_THROW(Dune::NotImplemented, "Internal energies!"); + DUNE_THROW(Dune::NotImplemented, "Enthalpy"); }; /*! @@ -319,7 +319,7 @@ public: int phaseIdx) { // TODO! - DUNE_THROW(Dune::NotImplemented, "Thermal conductivity!"); + DUNE_THROW(Dune::NotImplemented, "Thermal conductivity"); } /*! @@ -336,7 +336,7 @@ public: int phaseIdx) { // TODO! - DUNE_THROW(Dune::NotImplemented, "Heat capacity!"); + DUNE_THROW(Dune::NotImplemented, "Heat capacity"); } }; diff --git a/dumux/boxmodels/1p2c/1p2cvolumevariables.hh b/dumux/boxmodels/1p2c/1p2cvolumevariables.hh index 55da13154b622dfe10d320bf8ada621c83aa330a..4d41328539c3aa1ba778d1193651b78496c62e86 100644 --- a/dumux/boxmodels/1p2c/1p2cvolumevariables.hh +++ b/dumux/boxmodels/1p2c/1p2cvolumevariables.hh @@ -120,7 +120,6 @@ public: } fluidState_.setMoleFraction(/*phaseIdx=*/0, /*compIdx=*/0, 1 - x1); fluidState_.setMoleFraction(/*phaseIdx=*/0, /*compIdx=*/1, x1); - fluidState_.updateAverageMolarMass(/*phaseIdx=*/0); typename FluidSystem::ParameterCache paramCache; paramCache.updatePhase(fluidState_, /*phaseIdx=*/0); diff --git a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh index f7ab0f41d454fb71fa68107afe24379e3b103dcb..b3d040276c35b3be0cbeebb952c55d6d72d890d5 100644 --- a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh +++ b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh @@ -132,9 +132,9 @@ protected: { // copmute and set the internal energies of the fluid phases for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - Scalar u = FluidSystem::internalEnergy(this->fluidState_, paramCache, phaseIdx); + Scalar h = FluidSystem::enthalpy(this->fluidState_, paramCache, phaseIdx); - this->fluidState_.setInternalEnergy(phaseIdx, u); + this->fluidState_.setEnthalpy(phaseIdx, h); } // copmute and set the heat capacity of the solid phase diff --git a/dumux/boxmodels/2pni/2pnivolumevariables.hh b/dumux/boxmodels/2pni/2pnivolumevariables.hh index 297c6d53b824310f721ed18ec8f41e1227e1629d..5399ac08ac1a3a9abeedf1ab6843a2a700313e8e 100644 --- a/dumux/boxmodels/2pni/2pnivolumevariables.hh +++ b/dumux/boxmodels/2pni/2pnivolumevariables.hh @@ -129,9 +129,9 @@ protected: { // copmute and set the internal energies of the fluid phases for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - Scalar u = FluidSystem::internalEnergy(this->fluidState_, paramCache, phaseIdx); + Scalar h = FluidSystem::enthalpy(this->fluidState_, paramCache, phaseIdx); - this->fluidState_.setInternalEnergy(phaseIdx, u); + this->fluidState_.setEnthalpy(phaseIdx, h); } // copmute and set the heat capacity of the solid phase diff --git a/dumux/boxmodels/MpNc/MpNcvolumevariables.hh b/dumux/boxmodels/MpNc/MpNcvolumevariables.hh index 37038c8c4b1796431d0fcf4ef414b13614842950..163225029b4ceb11ebffc1f14f50cc434b0d014c 100644 --- a/dumux/boxmodels/MpNc/MpNcvolumevariables.hh +++ b/dumux/boxmodels/MpNc/MpNcvolumevariables.hh @@ -219,7 +219,6 @@ public: // update the remaining parts of the energy module EnergyVolumeVariables::update(fluidState_, paramCache, - priVars, element, elemGeom, scvIdx, diff --git a/dumux/boxmodels/MpNc/energy/MpNcfluxvariablesenergy.hh b/dumux/boxmodels/MpNc/energy/MpNcfluxvariablesenergy.hh index d4589e3b4c0899605f3717ccb5a793948a9befa0..82334764273967ea4b56b35519e7e690dfa0b71b 100644 --- a/dumux/boxmodels/MpNc/energy/MpNcfluxvariablesenergy.hh +++ b/dumux/boxmodels/MpNc/energy/MpNcfluxvariablesenergy.hh @@ -112,7 +112,7 @@ public: for (int scvIdx = 0; scvIdx < fvElemGeom.numVertices; scvIdx++) { tmp = fvElemGeom.subContVolFace[scvfIdx].grad[scvIdx]; - tmp *= elemVolVars[scvIdx].temperature(); + tmp *= elemVolVars[scvIdx].fluidState().temperature(/*phaseIdx=*/0); temperatureGradient += tmp; } @@ -121,14 +121,16 @@ public: lambdaPm_ = lumpedLambdaPm(problem, + element, fvElemGeom, scvfIdx, elemVolVars) ; } - Scalar lumpedLambdaPm(const Problem & problem, - const FVElementGeometry & fvElemGeom, + Scalar lumpedLambdaPm(const Problem &problem, + const Element &element, + const FVElementGeometry & fvElemGeom, int faceIdx, const ElementVolumeVariables & elemVolVars) { @@ -136,38 +138,46 @@ public: const int i = fvElemGeom.subContVolFace[faceIdx].i; const int j = fvElemGeom.subContVolFace[faceIdx].j; - const Scalar Sli = elemVolVars[i].fluidState().saturation(lPhaseIdx); - const Scalar Slj = elemVolVars[j].fluidState().saturation(lPhaseIdx); + typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables))::FluidState FluidState; + const FluidState &fsI = elemVolVars[i].fluidState(); + const FluidState &fsJ = elemVolVars[j].fluidState(); + const Scalar Sli = fsI.saturation(lPhaseIdx); + const Scalar Slj = fsJ.saturation(lPhaseIdx); - const Scalar Sl = std::max<Scalar>(0.0, 0.5*(Sli + Slj)); + typename FluidSystem::ParameterCache paramCacheI, paramCacheJ; + paramCacheI.updateAll(fsI); + paramCacheJ.updateAll(fsJ); - // const Scalar lambdaDry = 0.583; // W / (K m) // works, orig - // const Scalar lambdaWet = 1.13; // W / (K m) // works, orig - - const typename FluidSystem::MutableParameters mutParams; //dummy - const Scalar lambdaDry = 0.5 * (problem.spatialParameters().soilThermalConductivity() + FluidSystem::thermalConductivity(mutParams, gPhaseIdx) ); // W / (K m) - const Scalar lambdaWet = 0.5 * (problem.spatialParameters().soilThermalConductivity() + FluidSystem::thermalConductivity(mutParams, lPhaseIdx)) ; // W / (K m) - - // the heat conductivity of the matrix. in general this is a - // tensorial value, but we assume isotropic heat conductivity. - // This is the Sommerton approach with lambdaDry = - // lambdaSn100%. Taken from: H. Class: "Theorie und - // numerische Modellierung nichtisothermer Mehrphasenprozesse - // in NAPL-kontaminierten poroesen Medien", PhD Thesis, University of - // Stuttgart, Institute of Hydraulic Engineering, p. 57 - - Scalar result; - if (Sl < 0.1) { - // regularization - Dumux::Spline<Scalar> sp(0, 0.1, // x1, x2 - 0, sqrt(0.1), // y1, y2 - 5*0.5/sqrt(0.1), 0.5/sqrt(0.1)); // m1, m2 - result = lambdaDry + sp.eval(Sl)*(lambdaWet - lambdaDry); - } - else - result = lambdaDry + std::sqrt(Sl)*(lambdaWet - lambdaDry); + const Scalar Sl = std::max<Scalar>(0.0, 0.5*(Sli + Slj)); - return result; + // const Scalar lambdaDry = 0.583; // W / (K m) // works, orig + // const Scalar lambdaWet = 1.13; // W / (K m) // works, orig + + Scalar lambdaSoilI = problem.spatialParameters().soilThermalConductivity(element, fvElemGeom, i); + Scalar lambdaSoilJ = problem.spatialParameters().soilThermalConductivity(element, fvElemGeom, i); + const Scalar lambdaDry = 0.5 * (lambdaSoilI + FluidSystem::thermalConductivity(fsI, paramCacheI, gPhaseIdx)); // W / (K m) + const Scalar lambdaWet = 0.5 * (lambdaSoilJ + FluidSystem::thermalConductivity(fsJ, paramCacheJ, lPhaseIdx)) ; // W / (K m) + + // the heat conductivity of the matrix. in general this is a + // tensorial value, but we assume isotropic heat conductivity. + // This is the Sommerton approach with lambdaDry = + // lambdaSn100%. Taken from: H. Class: "Theorie und + // numerische Modellierung nichtisothermer Mehrphasenprozesse + // in NAPL-kontaminierten poroesen Medien", PhD Thesis, University of + // Stuttgart, Institute of Hydraulic Engineering, p. 57 + + Scalar result; + if (Sl < 0.1) { + // regularization + Dumux::Spline<Scalar> sp(0, 0.1, // x1, x2 + 0, sqrt(0.1), // y1, y2 + 5*0.5/sqrt(0.1), 0.5/sqrt(0.1)); // m1, m2 + result = lambdaDry + sp.eval(Sl)*(lambdaWet - lambdaDry); + } + else + result = lambdaDry + std::sqrt(Sl)*(lambdaWet - lambdaDry); + + return result; } /*! diff --git a/dumux/boxmodels/MpNc/energy/MpNclocalresidualenergy.hh b/dumux/boxmodels/MpNc/energy/MpNclocalresidualenergy.hh index cae84ef11739c749447285fcf17ff32bb32aac54..45e668dd49ba7bc95bca3f0a1c863cb73638926c 100644 --- a/dumux/boxmodels/MpNc/energy/MpNclocalresidualenergy.hh +++ b/dumux/boxmodels/MpNc/energy/MpNclocalresidualenergy.hh @@ -142,7 +142,7 @@ public: // heat stored in the rock matrix result[energyEqIdx] += - volVars.temperature() + volVars.fluidState().temperature(/*phaseIdx=*/0) * volVars.soilDensity() * (1.0 - volVars.porosity()) * volVars.heatCapacity(); diff --git a/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh b/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh index 7cd6b7681de56cb4508018111ea3e4ec3cf75991..e9bc865e628b0d945c630d411ee3c55fb2854967 100644 --- a/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh +++ b/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh @@ -88,7 +88,6 @@ public: */ void update(FluidState &fs, ParameterCache ¶mCache, - const PrimaryVariables &sol, const Element &element, const FVElementGeometry &elemGeom, int scvIdx, @@ -174,9 +173,8 @@ public: // set the enthalpies for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - Scalar u = - FluidSystem::internalEnergy(fs, paramCache, phaseIdx); - fs.setInternalEnergy(phaseIdx, u); + Scalar h = FluidSystem::enthalpy(fs, paramCache, phaseIdx); + fs.setEnthalpy(phaseIdx, h); } } diff --git a/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh b/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh index 594a9666797c2cc06e5176580b90984fa5cfad20..79cf8f60ed9d3bb44d35aef2704278c783d51166 100644 --- a/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh +++ b/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh @@ -26,7 +26,6 @@ #ifndef DUMUX_COMPOSITION_FROM_FUGACITIES_HH #define DUMUX_COMPOSITION_FROM_FUGACITIES_HH -#include "../MpNcfluidstates/nonequilibriumfluidstate.hh" #include <dumux/common/exceptions.hh> namespace Dumux { @@ -105,7 +104,6 @@ public: // right hand side Dune::FieldVector<Scalar, numComponents> b; - fluidState.updateAverageMolarMass(phaseIdx); paramCache.updatePhase(fluidState, phaseIdx); // maximum number of iterations @@ -195,7 +193,6 @@ protected: fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma); }; - fluidState.updateAverageMolarMass(phaseIdx); paramCache.updatePhase(fluidState, phaseIdx); Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx); @@ -240,7 +237,6 @@ protected: // deviate the mole fraction of the i-th component Scalar x_i = fluidState.moleFraction(phaseIdx, i); fluidState.setMoleFraction(phaseIdx, i, x_i + eps); - fluidState.updateAverageMolarMass(phaseIdx); paramCache.updatePhaseSingleMoleFraction(fluidState, phaseIdx, i); // compute new defect and derivative for all component @@ -266,7 +262,6 @@ protected: // reset composition to original value fluidState.setMoleFraction(phaseIdx, i, x_i); - fluidState.updateAverageMolarMass(phaseIdx); paramCache.updatePhaseSingleMoleFraction(fluidState, phaseIdx, i); // end forward differences @@ -325,7 +320,6 @@ protected: //sumMoleFrac += std::abs(newx); } - fluidState.updateAverageMolarMass(phaseIdx); paramCache.updatePhaseComposition(fluidState, phaseIdx); /* diff --git a/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh b/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh index 866b9d7c25b15faee377593cde1228b21408f837..f251dc06c051d2647c23abbbd87cb5a193f56894 100644 --- a/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh +++ b/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh @@ -30,7 +30,6 @@ #ifndef DUMUX_COMPUTE_FROM_REFERENCE_PHASE_HH #define DUMUX_COMPUTE_FROM_REFERENCE_PHASE_HH -#include "../MpNcfluidstates/genericfluidstate.hh" #include <dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh> namespace Dumux { @@ -107,24 +106,23 @@ public: ParameterCache ¶mCache, int refPhaseIdx, bool setViscosity, - bool setInternalEnergy) + bool setEnthalpy) { ComponentVector fugVec; // compute the density and enthalpy of the // reference phase - fluidState.updateAverageMolarMass(refPhaseIdx); paramCache.updatePhase(fluidState, refPhaseIdx); fluidState.setDensity(refPhaseIdx, FluidSystem::density(fluidState, paramCache, refPhaseIdx)); - if (setInternalEnergy) - fluidState.setInternalEnergy(refPhaseIdx, - FluidSystem::internalEnergy(fluidState, - paramCache, - refPhaseIdx)); + if (setEnthalpy) + fluidState.setEnthalpy(refPhaseIdx, + FluidSystem::enthalpy(fluidState, + paramCache, + refPhaseIdx)); if (setViscosity) fluidState.setViscosity(refPhaseIdx, @@ -157,11 +155,11 @@ public: paramCache, phaseIdx)); - if (setInternalEnergy) - fluidState.setInternalEnergy(phaseIdx, - FluidSystem::internalEnergy(fluidState, - paramCache, - phaseIdx)); + if (setEnthalpy) + fluidState.setEnthalpy(phaseIdx, + FluidSystem::enthalpy(fluidState, + paramCache, + phaseIdx)); } }; }; diff --git a/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh b/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh index 68b97fdf45d3c6174bf6769b79a6ae1bbcce83e8..d7db5a55124a62f691416ce4a119e022b84faabf 100644 --- a/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh +++ b/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh @@ -167,8 +167,6 @@ public: int rowIdx = phaseIdx*numComponents + compIdx; fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]); } - - fluidState.updateAverageMolarMass(phaseIdx); paramCache.updatePhaseComposition(fluidState, phaseIdx); Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx); @@ -180,8 +178,8 @@ public: } if (setInternalEnergy) { - value = FluidSystem::internalEnergy(fluidState, paramCache, phaseIdx); - fluidState.setInternalEnergy(phaseIdx, value); + value = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx); + fluidState.setEnthalpy(phaseIdx, value); } } }; diff --git a/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh b/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh index c44b79d1befbd3de13c10915f38b6c2bd838fce0..6568ec3317adc99f840fb17afadd0ed0956a388d 100644 --- a/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh +++ b/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh @@ -44,7 +44,19 @@ public: enum { numComponents = FluidSystem::numComponents }; EquilibriumFluidState() - { Valgrind::SetUndefined(*this); } + { + // set the composition to 0 + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + for (int compIdx = 0; compIdx < numComponents; ++compIdx) + moleFraction_[phaseIdx][compIdx] = 0; + + averageMolarMass_[phaseIdx] = 0; + sumMoleFractions_[phaseIdx] = 0; + } + + // make everything undefined so that valgrind will complain + Valgrind::SetUndefined(*this); + } template <class FluidState> EquilibriumFluidState(FluidState &fs) @@ -72,29 +84,12 @@ public: Scalar massFraction(int phaseIdx, int compIdx) const { return - sumMassFractions(phaseIdx)* - molarity(phaseIdx, compIdx)* - FluidSystem::molarMass(compIdx) - / - density(phaseIdx); + sumMoleFractions_[phaseIdx] + * moleFraction_[phaseIdx][compIdx] + * FluidSystem::molarMass(compIdx) + / averageMolarMass_[phaseIdx]; } - /*! - * \brief The sum of all component mole fractions in a phase [] - * - * We define this to be the same as the sum of all mass fractions. - */ - Scalar sumMoleFractions(int phaseIdx) const - { return sumMoleFractions_[phaseIdx]; } - - /*! - * \brief The sum of all component mass fractions in a phase [] - * - * We define this to be the same as the sum of all mole fractions. - */ - Scalar sumMassFractions(int phaseIdx) const - { return sumMoleFractions_[phaseIdx]; } - /*! * \brief The average molar mass of a fluid phase [kg/mol] */ @@ -159,13 +154,13 @@ public: * \brief The specific enthalpy of a fluid phase [J/kg] */ Scalar enthalpy(int phaseIdx) const - { return internalEnergy_[phaseIdx] + pressure(phaseIdx)/(density(phaseIdx)); } + { return enthalpy_[phaseIdx]; } /*! * \brief The specific internal energy of a fluid phase [J/kg] */ Scalar internalEnergy(int phaseIdx) const - { return internalEnergy_[phaseIdx]; } + { return enthalpy_[phaseIdx] - pressure(phaseIdx)/density(phaseIdx); } /*! * \brief The dynamic viscosity of a fluid phase [Pa s] @@ -185,16 +180,6 @@ public: Scalar temperature() const { return temperature_; } - /*! - * \brief The capillary pressure [Pa] between a phase and a - * reference phase. - * - * In order to make the term "capillary pressure" meaningful in a - * physical sense, mechanic equilibrium needs to be assumed. - */ - Scalar capillaryPressure(int refPhaseIdx, int phaseIdx) const - { return pressure_[phaseIdx] - pressure_[refPhaseIdx]; } - /*! * \brief The fugacity of a component * @@ -230,7 +215,7 @@ public: pressure_[phaseIdx] = fs.pressure(phaseIdx); saturation_[phaseIdx] = fs.saturation(phaseIdx); density_[phaseIdx] = fs.density(phaseIdx); - internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx); + enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx); viscosity_[phaseIdx] = fs.viscosity(phaseIdx); } temperature_ = fs.temperature(0); @@ -258,7 +243,19 @@ public: * \brief Set the mole fraction of a component in a phase [] */ void setMoleFraction(int phaseIdx, int compIdx, Scalar value) - { moleFraction_[phaseIdx][compIdx] = value; } + { + Valgrind::CheckDefined(value); + + Scalar delta = value - moleFraction_[phaseIdx][compIdx]; + + moleFraction_[phaseIdx][compIdx] = value; + + sumMoleFractions_[phaseIdx] += delta; + averageMolarMass_[phaseIdx] += delta*FluidSystem::molarMass(compIdx); + + Valgrind::SetDefined(sumMoleFractions_[phaseIdx]); + Valgrind::SetDefined(averageMolarMass_[phaseIdx]); + }; /*! * \brief Set the fugacity of a component in a phase [] @@ -275,8 +272,8 @@ public: /*! * \brief Set the specific enthalpy of a phase [J/m^3] */ - void setInternalEnergy(int phaseIdx, Scalar value) - { internalEnergy_[phaseIdx] = value; } + void setEnthalpy(int phaseIdx, Scalar value) + { enthalpy_[phaseIdx] = value; } /*! * \brief Set the dynamic viscosity of a phase [Pa s] @@ -284,24 +281,6 @@ public: void setViscosity(int phaseIdx, Scalar value) { viscosity_[phaseIdx] = value; } - /*! - * \brief Calculatate the mean molar mass of a phase given that - * all mole fractions have been set - */ - void updateAverageMolarMass(int phaseIdx) - { - averageMolarMass_[phaseIdx] = 0; - sumMoleFractions_[phaseIdx] = 0; - - for (int compIdx = 0; compIdx < numComponents; ++compIdx) { - sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx]; - averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx); - } - sumMoleFractions_[phaseIdx] = std::max(1e-15, sumMoleFractions_[phaseIdx]); - Valgrind::CheckDefined(averageMolarMass_[phaseIdx]); - Valgrind::CheckDefined(sumMoleFractions_[phaseIdx]); - } - /*! * \brief Make sure that all attributes are defined. * @@ -339,7 +318,7 @@ protected: Scalar pressure_[numPhases]; Scalar saturation_[numPhases]; Scalar density_[numPhases]; - Scalar internalEnergy_[numPhases]; + Scalar enthalpy_[numPhases]; Scalar viscosity_[numPhases]; Scalar temperature_; }; diff --git a/dumux/material/MpNcfluidstates/immisciblefluidstate.hh b/dumux/material/MpNcfluidstates/immisciblefluidstate.hh index 69836d04a53cdd07f70ff77e68962e0e8a31d97e..eb087e6e70d0acba90521754fa8977ce58fe8b69 100644 --- a/dumux/material/MpNcfluidstates/immisciblefluidstate.hh +++ b/dumux/material/MpNcfluidstates/immisciblefluidstate.hh @@ -74,22 +74,6 @@ public: Scalar massFraction(int phaseIdx, int compIdx) const { return (phaseIdx == compIdx)?1.0:0.0; } - /*! - * \brief The sum of all component mole fractions in a phase [] - * - * We define this to be the same as the sum of all mass fractions. - */ - Scalar sumMoleFractions(int phaseIdx) const - { return 1.0; } - - /*! - * \brief The sum of all component mass fractions in a phase [] - * - * We define this to be the same as the sum of all mole fractions. - */ - Scalar sumMassFractions(int phaseIdx) const - { return 1.0; } - /*! * \brief The average molar mass of a fluid phase [kg/mol] */ @@ -175,13 +159,13 @@ public: * \brief The specific enthalpy of a fluid phase [J/kg] */ Scalar enthalpy(int phaseIdx) const - { return internalEnergy_[phaseIdx] + pressure(phaseIdx)/(density(phaseIdx)); } + { return enthalpy_[phaseIdx]; } /*! * \brief The specific internal energy of a fluid phase [J/kg] */ Scalar internalEnergy(int phaseIdx) const - { return internalEnergy_[phaseIdx]; } + { return enthalpy_[phaseIdx] - pressure(phaseIdx)/density(phaseIdx); } /*! * \brief The dynamic viscosity of a fluid phase [Pa s] @@ -189,7 +173,6 @@ public: Scalar viscosity(int phaseIdx) const { return viscosity_[phaseIdx]; } - /***************************************************** * Access to fluid properties which only make sense * if assuming thermodynamic equilibrium @@ -231,7 +214,7 @@ public: pressure_[phaseIdx] = fs.pressure(phaseIdx); saturation_[phaseIdx] = fs.saturation(phaseIdx); density_[phaseIdx] = fs.density(phaseIdx); - internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx); + enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx); viscosity_[phaseIdx] = fs.viscosity(phaseIdx); } temperature_ = fs.temperature(0); @@ -262,10 +245,10 @@ public: { density_[phaseIdx] = value; } /*! - * \brief Set the specific internal energy of a phase [J/m^3] + * \brief Set the specific enthalpy of a phase [J/m^3] */ - void setInternalEnergy(int phaseIdx, Scalar value) - { internalEnergy_[phaseIdx] = value; } + void setEnthalpy(int phaseIdx, Scalar value) + { enthalpy_[phaseIdx] = value; } /*! * \brief Set the dynamic viscosity of a phase [Pa s] @@ -288,7 +271,7 @@ public: //for (int j = 0; j < numComponents; ++j) { // Valgrind::CheckDefined(fugacityCoefficient_[i][j]); //} - Valgrind::CheckDefined(pressure_[i]); + Valgrind:CheckDefined(pressure_[i]); Valgrind::CheckDefined(saturation_[i]); Valgrind::CheckDefined(density_[i]); //Valgrind::CheckDefined(internalEnergy_[i]); @@ -303,7 +286,7 @@ protected: Scalar pressure_[numPhases]; Scalar saturation_[numPhases]; Scalar density_[numPhases]; - Scalar internalEnergy_[numPhases]; + Scalar enthalpy_[numPhases]; Scalar viscosity_[numPhases]; Scalar temperature_; }; diff --git a/dumux/material/MpNcfluidstates/genericfluidstate.hh b/dumux/material/MpNcfluidstates/nonequilibriumfluidstate.hh similarity index 83% rename from dumux/material/MpNcfluidstates/genericfluidstate.hh rename to dumux/material/MpNcfluidstates/nonequilibriumfluidstate.hh index 9b4778c7c3033a9244627842cdd637ba3fdb1777..c6e64b4bdf6c844a471b297af12e299be7943e29 100644 --- a/dumux/material/MpNcfluidstates/genericfluidstate.hh +++ b/dumux/material/MpNcfluidstates/nonequilibriumfluidstate.hh @@ -40,14 +40,26 @@ namespace Dumux * any assumptions. */ template <class Scalar, class FluidSystem> -class GenericFluidState +class NonEquilibriumFluidState { public: enum { numPhases = FluidSystem::numPhases }; enum { numComponents = FluidSystem::numComponents }; - GenericFluidState() - { Valgrind::SetUndefined(*this); } + NonEquilibriumFluidState() + { + // set the composition to 0 + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + for (int compIdx = 0; compIdx < numComponents; ++compIdx) + moleFractions_[phaseIdx][compIdx] = 0; + + averageMolarMass_[phaseIdx] = 0; + sumMoleFractions_[phaseIdx] = 0; + } + + // make everything undefined so that valgrind will complain + Valgrind::SetUndefined(*this); + } /***************************************************** * Generic access to fluid properties (No assumptions @@ -72,28 +84,12 @@ public: Scalar massFraction(int phaseIdx, int compIdx) const { return - sumMassFraction(phaseIdx) + sumMoleFractions_[phaseIdx] * moleFraction_[phaseIdx][compIdx] * FluidSystem::molarMass(compIdx) / averageMolarMass_[phaseIdx]; } - /*! - * \brief The sum of all component mole fractions in a phase [] - * - * We define this to be the same as the sum of all mass fractions. - */ - Scalar sumMoleFractions(int phaseIdx) const - { return sumMoleFractions_[phaseIdx]; } - - /*! - * \brief The sum of all component mass fractions in a phase [] - * - * We define this to be the same as the sum of all mole fractions. - */ - Scalar sumMassFraction(int phaseIdx) const - { return sumMoleFractions_[phaseIdx]; } - /*! * \brief The mean molar mass of a fluid phase [kg/mol] */ @@ -153,18 +149,17 @@ public: */ Scalar pressure(int phaseIdx) const { return pressure_[phaseIdx]; }; - /*! - * \brief The specific internal energy of a fluid phase [J/kg] + * \brief The specific enthalpy of a fluid phase [J/kg] */ - Scalar internalEnergy(int phaseIdx) const - { return internalEnergy_[phaseIdx]; }; + Scalar enthalpy(int phaseIdx) const + { return enthalpy_[phaseIdx]; } /*! - * \brief The specific enthalpy of a fluid phase [J/kg] + * \brief The specific internal energy of a fluid phase [J/kg] */ - Scalar enthalpy(int phaseIdx) const - { return internalEnergy(phaseIdx) + pressure(phaseIdx)/density(phaseIdx); }; + Scalar internalEnergy(int phaseIdx) const + { return enthalpy_[phaseIdx] - pressure(phaseIdx)/density(phaseIdx); } /*! * \brief The dynamic viscosity of a fluid phase [Pa s] @@ -199,7 +194,19 @@ public: * \brief Set the mole fraction of a component in a phase [] */ void setMoleFraction(int phaseIdx, int compIdx, Scalar value) - { moleFraction_[phaseIdx][compIdx] = value; } + { + Valgrind::CheckDefined(value); + + Scalar delta = value - moleFraction_[phaseIdx][compIdx]; + + moleFraction_[phaseIdx][compIdx] = value; + + sumMoleFractions_[phaseIdx] += delta; + averageMolarMass_[phaseIdx] += delta*FluidSystem::molarMass(compIdx); + + Valgrind::SetDefined(sumMolarMass_[phaseIdx]); + Valgrind::SetDefined(averageMolarMass_[phaseIdx]); + } /*! * \brief Set the fugacity of a component in a phase [] @@ -214,10 +221,10 @@ public: { density_[phaseIdx] = value; } /*! - * \brief Set the specific internal energy of a phase [J/m^3] + * \brief Set the specific enthalpy of a phase [J/m^3] */ - void setInternalEnergy(int phaseIdx, Scalar value) - { internalEnergy_[phaseIdx] = value; } + void setEnthalpy(int phaseIdx, Scalar value) + { enthalpy_[phaseIdx] = value; } /*! * \brief Set the dynamic viscosity of a phase [Pa s] @@ -225,24 +232,6 @@ public: void setViscosity(int phaseIdx, Scalar value) { viscosity_[phaseIdx] = value; } - /*! - * \brief Calculatate the mean molar mass of a phase given that - * all mole fractions have been set - */ - void updateAverageMolarMass(int phaseIdx) - { - averageMolarMass_[phaseIdx] = 0; - sumMoleFractions_[phaseIdx] = 0; - - for (int compIdx = 0; compIdx < numComponents; ++compIdx) { - sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx]; - averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx); - } - sumMoleFractions_[phaseIdx] = std::max(1e-15, sumMoleFractions_[phaseIdx]); - Valgrind::CheckDefined(averageMolarMass_[phaseIdx]); - Valgrind::CheckDefined(sumMoleFractions_[phaseIdx]); - } - /*! * \brief Retrieve all parameters from an arbitrary fluid * state. @@ -260,7 +249,7 @@ public: pressure_[phaseIdx] = fs.pressure(phaseIdx); saturation_[phaseIdx] = fs.saturation(phaseIdx); density_[phaseIdx] = fs.density(phaseIdx); - internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx); + enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx); viscosity_[phaseIdx] = fs.viscosity(phaseIdx); temperature_[phaseIdx] = fs.temperature(phaseIdx); } @@ -302,7 +291,7 @@ protected: Scalar pressure_[numPhases]; Scalar saturation_[numPhases]; Scalar density_[numPhases]; - Scalar internalEnergy_[numPhases]; + Scalar enthalpy_[numPhases]; Scalar viscosity_[numPhases]; Scalar temperature_[numPhases]; }; diff --git a/dumux/material/MpNcfluidsystems/1pfluidsystem.hh b/dumux/material/MpNcfluidsystems/1pfluidsystem.hh index bac653903f6bb7ba21c57855234ff290b2642c44..9fcd1a85ec911ae36e5e24a2090a87c30763c47b 100644 --- a/dumux/material/MpNcfluidsystems/1pfluidsystem.hh +++ b/dumux/material/MpNcfluidsystems/1pfluidsystem.hh @@ -303,22 +303,22 @@ public: /*! * \brief Given a phase's composition, temperature, pressure and - * density, calculate its specific internal energy [J/kg]. + * density, calculate its specific enthalpy [J/kg]. * * \param fluidState An abitrary fluid state * \param paramCache The fluid system's parameter cache * \param phaseIdx The index of the fluid phase to consider */ template <class FluidState> - static Scalar internalEnergy(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + static Scalar enthalpy(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { assert(0 <= phaseIdx && phaseIdx < numPhases); Scalar temperature = fluidState.temperature(phaseIdx); Scalar pressure = fluidState.pressure(phaseIdx); - return Fluid::internalEnergy(temperature, pressure); + return Fluid::enthalpy(temperature, pressure); } /*! diff --git a/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh b/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh index dc36667a279578a7a320ea98791de16b798437c0..70196ee0e1a44f419eb3eca93fcc0afec0209257 100644 --- a/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh +++ b/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh @@ -329,18 +329,11 @@ public: }; /*! - * \brief Return the specific internal of a fluid phase [J/kg]. - * - * \param phaseIdx index of the phase - * \param temperature phase temperature in [K] - * \param pressure phase pressure in [Pa] - * \param fluidState The fluid state of the two-phase model - * \tparam FluidState the fluid state class of the two-phase model - * \return returns the specific enthalpy of the phase [J/kg] + * \brief Return the specific enthalpy of a fluid phase [J/kg]. */ using Base::internalEnergy; template <class FluidState> - static Scalar internalEnergy(const FluidState &fluidState, + static Scalar enthalpy(const FluidState &fluidState, int phaseIdx) { assert(0 <= phaseIdx && phaseIdx < numPhases); @@ -348,8 +341,8 @@ public: Scalar temperature = fluidState.temperature(phaseIdx); Scalar pressure = fluidState.pressure(phaseIdx); if (phaseIdx == wPhaseIdx) - return WettingPhase::internalEnergy(temperature, pressure); - return NonWettingPhase::internalEnergy(temperature, pressure); + return WettingPhase::enthalpy(temperature, pressure); + return NonWettingPhase::enthalpy(temperature, pressure); } /*! diff --git a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh index f643546ed46699cf8783c25d43e97ccb4c7febc0..e5b4f0ed388e51183efbc499aa9e0a163d3c659a 100644 --- a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh +++ b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh @@ -25,8 +25,6 @@ #ifndef DUMUX_H2O_N2_FLUID_SYSTEM_HH #define DUMUX_H2O_N2_FLUID_SYSTEM_HH -#include <dumux/material/MpNcfluidstates/nonequilibriumfluidstate.hh> - #include <dumux/material/components/simpleh2o.hh> #include <dumux/material/components/h2o.hh> #include <dumux/material/components/n2.hh> @@ -34,7 +32,6 @@ #include <dumux/material/idealgas.hh> #include <dumux/material/binarycoefficients/h2o_n2.hh> -#include <dumux/material/MpNcfluidstates/nonequilibriumfluidstate.hh> #include <dumux/common/valgrind.hh> #include <dumux/common/exceptions.hh> @@ -456,7 +453,7 @@ public: /*! * \brief Given a phase's composition, temperature, pressure and - * density, calculate its specific internal energy [J/kg]. + * density, calculate its specific enthalpy [J/kg]. * * \todo This fluid system neglects the contribution of * gas-molecules in the liquid phase. This contribution is @@ -468,9 +465,9 @@ public: * \param phaseIdx The index of the fluid phase to consider */ template <class FluidState> - static Scalar internalEnergy(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + static Scalar enthalpy(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { Scalar T = fluidState.temperature(phaseIdx); Scalar p = fluidState.pressure(phaseIdx); @@ -479,15 +476,15 @@ public: if (phaseIdx == lPhaseIdx) { // TODO: correct way to deal with the solutes??? return - H2O::liquidInternalEnergy(T, p) ; + H2O::liquidEnthalpy(T, p) ; } else { // assume ideal gas Scalar XH2O = fluidState.massFraction(gPhaseIdx, H2OIdx); Scalar XN2 = fluidState.massFraction(gPhaseIdx, N2Idx); Scalar result = 0; - result += XH2O*H2O::gasInternalEnergy(T, p); - result += XN2*N2::gasInternalEnergy(T, p); + result += XH2O*H2O::gasEnthalpy(T, p); + result += XN2*N2::gasEnthalpy(T, p); return result; } } diff --git a/test/boxmodels/1p2c/interstitialfluidtrailfluidsystem.hh b/test/boxmodels/1p2c/interstitialfluidtrailfluidsystem.hh index a7aa775159520e8cd3103e1a3a69b1522a246013..ab2b6bd91e581a96cc594ca84d189b2af67903a8 100644 --- a/test/boxmodels/1p2c/interstitialfluidtrailfluidsystem.hh +++ b/test/boxmodels/1p2c/interstitialfluidtrailfluidsystem.hh @@ -301,13 +301,13 @@ public: * \param phaseIdx for which phase to give back the heat capacity */ template <class FluidState> - static Scalar internalEnergy(const FluidState &fluidState, + static Scalar enthalpy(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx) { assert(0 <= phaseIdx && phaseIdx < numPhases); - DUNE_THROW(Dune::NotImplemented, "Internal energies"); + DUNE_THROW(Dune::NotImplemented, "Enthalpies"); }; /*! diff --git a/test/boxmodels/MpNc/obstaclespatialparameters.hh b/test/boxmodels/MpNc/obstaclespatialparameters.hh index 60a2d117030eb025437d65f2eab40b2524324983..ddbac95ebde0f4dc0b98ebd2663685d3c9daf95f 100644 --- a/test/boxmodels/MpNc/obstaclespatialparameters.hh +++ b/test/boxmodels/MpNc/obstaclespatialparameters.hh @@ -272,16 +272,16 @@ public: return coarseMaterialParams_; } - const Scalar & soilDensity(const Element &element, - const FVElementGeometry &fvElemGeom, - int scvIdx) const + Scalar soilDensity(const Element &element, + const FVElementGeometry &fvElemGeom, + int scvIdx) const { return 2700. ; // density of granite [kg/m^3] } - const Scalar & soilThermalConductivity(const Element &element, - const FVElementGeometry &fvElemGeom, - int scvIdx) const + Scalar soilThermalConductivity(const Element &element, + const FVElementGeometry &fvElemGeom, + int scvIdx) const { return 2.8; // conductivity of granite [W / (m K ) ] }