diff --git a/dumux/discretization/box/fourierslawnonequilibrium.hh b/dumux/discretization/box/fourierslawnonequilibrium.hh
index 5cb128deee82ab0ecf6b8109230d51780d0d0291..fb569e2a09a45721a792735296455cfb815cb01a 100644
--- a/dumux/discretization/box/fourierslawnonequilibrium.hh
+++ b/dumux/discretization/box/fourierslawnonequilibrium.hh
@@ -103,8 +103,8 @@ public:
//solid phase
else
{
- insideLambda += insideVolVars.solidThermalConductivity()*(1-insideVolVars.porosity());
- outsideLambda +=outsideVolVars.solidThermalConductivity()*(1-outsideVolVars.porosity());
+ insideLambda += insideVolVars.solidThermalConductivity()*(1.0-insideVolVars.porosity());
+ outsideLambda += outsideVolVars.solidThermalConductivity()*(1.0-outsideVolVars.porosity());
}
// scale by extrusion factor
diff --git a/dumux/discretization/cellcentered/tpfa/fourierslawnonequilibrium.hh b/dumux/discretization/cellcentered/tpfa/fourierslawnonequilibrium.hh
index f2c6e4acdbc67e21cb0e3bdd71435405919f2a7b..df203a6d334bfb11948ff7d4d1c6fab1cce2b808 100644
--- a/dumux/discretization/cellcentered/tpfa/fourierslawnonequilibrium.hh
+++ b/dumux/discretization/cellcentered/tpfa/fourierslawnonequilibrium.hh
@@ -95,7 +95,7 @@ public:
else //temp solid
{
tInside += elemVolVars[scvf.insideScvIdx()].temperatureSolid();
- tOutside += elemVolVars[scvf.outsideScvIdx()].temperatureSolid();
+ tOutside += elemVolVars[scvf.outsideScvIdx()].temperatureSolid();
}
Scalar tij = calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx);
@@ -134,7 +134,7 @@ public:
//solid phase
else
{
- insideLambda += insideVolVars.solidThermalConductivity()*(1-insideVolVars.porosity());
+ insideLambda += insideVolVars.solidThermalConductivity()*(1.0-insideVolVars.porosity());
}
const Scalar ti = computeTpfaTransmissibility(scvf, insideScv, insideLambda, insideVolVars.extrusionFactor());
@@ -167,7 +167,7 @@ public:
//solid phase
else
{
- outsideLambda +=outsideVolVars.solidThermalConductivity()*(1-outsideVolVars.porosity());
+ outsideLambda +=outsideVolVars.solidThermalConductivity()*(1.0-outsideVolVars.porosity());
}
Scalar tj;
if (dim == dimWorld)
diff --git a/dumux/material/fluidmatrixinteractions/1pia/fluidsolidinterfacialareashiwang.hh b/dumux/material/fluidmatrixinteractions/1pia/fluidsolidinterfacialareashiwang.hh
index ece673473bb93dd149988325de1c16bc8a217c02..77a9e6332dcbc683f17a0cf6242369de18c4bb71 100644
--- a/dumux/material/fluidmatrixinteractions/1pia/fluidsolidinterfacialareashiwang.hh
+++ b/dumux/material/fluidmatrixinteractions/1pia/fluidsolidinterfacialareashiwang.hh
@@ -45,7 +45,7 @@ public:
*/
static Scalar fluidSolidInterfacialArea(const Scalar porosity,
const Scalar characteristicLength)
- { return 6*(1-porosity)/characteristicLength; }
+ { return 6.0*(1.0-porosity)/characteristicLength; }
};
}
#endif
diff --git a/dumux/material/fluidmatrixinteractions/mp/mpadapter.hh b/dumux/material/fluidmatrixinteractions/mp/mpadapter.hh
index 698615cb5ea4386c568037b92b63873d765393a8..2754859c6d4e0463f783a6beeda8af106f22b0f3 100644
--- a/dumux/material/fluidmatrixinteractions/mp/mpadapter.hh
+++ b/dumux/material/fluidmatrixinteractions/mp/mpadapter.hh
@@ -28,6 +28,8 @@
#define DUMUX_MP_2P_ADAPTER_HH
#include <algorithm>
+#include <cassert>
+#include <dumux/common/typetraits/typetraits.hh>
namespace Dumux
{
@@ -58,10 +60,10 @@ public:
const FluidState &state,
const int wPhaseIdx)
{
+ assert(values.size() == 2);
const int nPhaseIdx = 1 - wPhaseIdx;
// non-wetting phase gets the capillary pressure added
values[nPhaseIdx] = 0;
-
// wetting phase does not get anything added
values[wPhaseIdx] = - MaterialLaw::pc(params, state.saturation(wPhaseIdx));
}
@@ -78,6 +80,7 @@ public:
const FluidState &state,
const int wPhaseIdx)
{
+ assert(values.size() == 2);
const int nPhaseIdx = 1 - wPhaseIdx;
values[wPhaseIdx] = MaterialLaw::krw(params, state.saturation(wPhaseIdx));
values[nPhaseIdx] = MaterialLaw::krn(params, state.saturation(wPhaseIdx));
diff --git a/dumux/material/spatialparams/fvnonequilibrium.hh b/dumux/material/spatialparams/fvnonequilibrium.hh
index 73ed50d919d37f82e38692b2b902ed7129fabf15..c2cf3be5563f787914c8d81f34fac0dbd8b6b443 100644
--- a/dumux/material/spatialparams/fvnonequilibrium.hh
+++ b/dumux/material/spatialparams/fvnonequilibrium.hh
@@ -34,17 +34,13 @@ namespace Dumux {
*/
template<class FVGridGeometry, class Scalar, class Implementation>
class FVNonEquilibriumSpatialParams
-: public FVSpatialParams<FVGridGeometry,
- Scalar,
- Implementation>
+: public FVSpatialParams<FVGridGeometry, Scalar, Implementation>
{
-
using ParentType = FVSpatialParams<FVGridGeometry, Scalar, Implementation>;
using GridView = typename FVGridGeometry::GridView;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVGridGeometry::SubControlVolume;
using Element = typename GridView::template Codim<0>::Entity;
-
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
diff --git a/dumux/porousmediumflow/1pnc/model.hh b/dumux/porousmediumflow/1pnc/model.hh
index 73c1b5adf83551de60c5b2e26e18d6769f9512a0..856bfdd71674bc125fa4fa379553fc69def5e5cd 100644
--- a/dumux/porousmediumflow/1pnc/model.hh
+++ b/dumux/porousmediumflow/1pnc/model.hh
@@ -249,7 +249,7 @@ public:
} // end namespace Properties
template<class OnePNCModelTraits>
-struct OnePNCNonequilibriumModelTraits : public OnePNCModelTraits
+struct OnePNCUnconstrainedModelTraits : public OnePNCModelTraits
{
static constexpr int numConstraintEq() { return 0; }
};
@@ -262,6 +262,7 @@ namespace TTag {
struct OnePNCNonEquil { using InheritsFrom = std::tuple<NonEquilibrium, OnePNC>; };
} // end namespace TTag
+
/////////////////////////////////////////////////
// Properties for the non-equilibrium OnePNC model
/////////////////////////////////////////////////
@@ -289,6 +290,7 @@ private:
public:
using type = NonEquilTraits;
};
+
// by default chemical non equilibrium is enabled in the nonequil model, switch that off here
template<class TypeTag>
struct EnableChemicalNonEquilibrium<TypeTag, TTag::OnePNCNonEquil> { static constexpr bool value = false; };
@@ -301,7 +303,7 @@ private:
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using EquilibriumTraits = OnePNCModelTraits<FluidSystem::numComponents, getPropValue<TypeTag, Properties::UseMoles>(), getPropValue<TypeTag, Properties::ReplaceCompEqIdx>()>;
public:
- using type = OnePNCNonequilibriumModelTraits<EquilibriumTraits>;
+ using type = OnePNCUnconstrainedModelTraits<EquilibriumTraits>;
};
//! in case we do not assume full non-equilibrium one needs a thermal conductivity
diff --git a/dumux/porousmediumflow/mpnc/volumevariables.hh b/dumux/porousmediumflow/mpnc/volumevariables.hh
index d2291213d8d0a3ddb3bb58b933c5c319f9ae5f35..1abb9cf00a012653375c3d26a7196fedaf8af436 100644
--- a/dumux/porousmediumflow/mpnc/volumevariables.hh
+++ b/dumux/porousmediumflow/mpnc/volumevariables.hh
@@ -188,7 +188,7 @@ public:
const auto& materialParams =
problem.spatialParams().materialLawParams(element, scv, elemSol);
// capillary pressures
- Scalar capPress[numPhases()];
+ std::vector<Scalar> capPress(numPhases());
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
using MPAdapter = MPAdapter<MaterialLaw, numPhases()>;
MPAdapter::capillaryPressures(capPress, materialParams, fluidState, wPhaseIdx);
@@ -494,7 +494,7 @@ protected:
std::array<std::array<Scalar, numFluidComps-1>, numPhases()> diffCoefficient_;
Scalar porosity_; //!< Effective porosity within the control volume
- Scalar relativePermeability_[numPhases()]; //!< Effective relative permeability within the control volume
+ std::array<Scalar, ModelTraits::numPhases()> relativePermeability_; //!< Effective relative permeability within the control volume
PermeabilityType permeability_;
//! Mass fractions of each component within each phase
@@ -568,9 +568,9 @@ public:
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
using MPAdapter = MPAdapter<MaterialLaw, numPhases()>;
MPAdapter::relativePermeabilities(relativePermeability_,
- materialParams,
- fluidState_,
- wPhaseIdx);
+ materialParams,
+ fluidState_,
+ wPhaseIdx);
typename FluidSystem::ParameterCache paramCache;
paramCache.updateAll(fluidState_);
if (enableDiffusion)
@@ -643,7 +643,7 @@ public:
problem.spatialParams().materialLawParams(element, scv, elemSol);
const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
// capillary pressures
- Scalar capPress[numPhases()];
+ std::vector<Scalar> capPress(numPhases());
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
using MPAdapter = MPAdapter<MaterialLaw, numPhases()>;
MPAdapter::capillaryPressures(capPress, materialParams, fluidState, wPhaseIdx);
@@ -1007,7 +1007,7 @@ protected:
}
std::array<std::array<Scalar, numFluidComps-1>, numPhases()> diffCoefficient_;
- Scalar relativePermeability_[numPhases()]; //!< Effective relative permeability within the control volume
+ std::array<Scalar, ModelTraits::numPhases()> relativePermeability_; //!< Effective relative permeability within the control volume
PermeabilityType permeability_;
Scalar xEquil_[numPhases()][numFluidComps];
diff --git a/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh b/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
index 01db8cbafb2827afdd4cacc76bbfa0a7ae824a28..b55b9e0e0ded869fc52c7c9b1c145237742e9531 100644
--- a/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
+++ b/dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh
@@ -148,9 +148,7 @@ public:
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);
- }
}
/*!
diff --git a/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh b/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
index b934ea5dd2f6040504b504b65201487e04937031..cfdbbda08ff585ff693cf871222811a2eb0e653e 100644
--- a/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
+++ b/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
@@ -456,7 +456,7 @@ private:
const auto &materialParams =
this->spatialParams().materialLawParamsAtPos(globalPos);
- Scalar capPress[numPhases];
+ std::vector<Scalar> capPress(numPhases);
//get the index for the wettingphase
const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
//obtain pc according to saturation
diff --git a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh
index 418c5dd4708fbb5c339d3df7944c9b7c038f61df..4810977f3f39625fe36d3318a3685eef3309e0d6 100644
--- a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh
+++ b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh
@@ -87,23 +87,22 @@ public:
const Scalar as = volVars.fluidSolidInterfacialArea();
//temperature fluid is the same for both fluids
- const Scalar TFluid = volVars.temperatureFluid(0);
- const Scalar TSolid = volVars.temperatureSolid();
+ 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 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");
+ 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() ;
@@ -112,7 +111,7 @@ public:
if (satW < (0 + eps) )
{
- solidToFluidEnergyExchange *= volVars.nusseltNumber(1) ;
+ solidToFluidEnergyExchange *= volVars.nusseltNumber(1) ;
}
else if ( (satW >= 0 + eps) and (satW < 1.0-eps) )
{
@@ -125,11 +124,11 @@ public:
QBoilFunc(volVars, 0.0), QBoilFunc(volVars, epsRegul), // y1, y2
0.0,dQBoil_dSw(volVars, epsRegul)); // m1, m2
- qBoil = sp.eval(satW) ;
+ qBoil = sp.eval(satW);
}
else if (satW>= (1.0-epsRegul) )
- {// regularize
+ { // regularize
typedef Dumux::Spline<Scalar> Spline;
Spline sp(1.0-epsRegul, 1.0, // x1, x2
QBoilFunc(volVars, 1.0-epsRegul), 0.0, // y1, y2
@@ -140,7 +139,7 @@ public:
else
{
qBoil = QBoilFunc(volVars, satW) ;
- }
+ }
solidToFluidEnergyExchange += qBoil;
}
@@ -149,8 +148,7 @@ public:
solidToFluidEnergyExchange *= volVars.nusseltNumber(0) ;
}
else
- DUNE_THROW(Dune::NotImplemented,
- "wrong range");
+ DUNE_THROW(Dune::NotImplemented, "wrong range");
using std::isfinite;
if (!isfinite(solidToFluidEnergyExchange))
@@ -184,27 +182,26 @@ public:
{
// 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() ;
+ const auto& fs = volVars.fluidState();
+ const Scalar g( 9.81 );
+ const Scalar gamma(0.0589);
+ const Scalar TSolid = volVars.temperatureSolid();
using std::pow;
- const Scalar as = 6.0 * (1.0-volVars.porosity()) / characteristicLength ;
- const Scalar mul = fs.viscosity(0) ;
+ const Scalar as = volVars.fluidSolidInterfacialArea();
+ 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 deltaRho = fs.density(0) - fs.density(1);
const Scalar firstBracket = pow(g * deltaRho / gamma, 0.5);
- const Scalar cp = FluidSystem::heatCapacity(fs, 0) ;
+ 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;
+ 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.
diff --git a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh
index 302a7d8f560fd79603f574553207b183adc64344..5784778fdaf8042dd73005ceef2d5acbe6a3a023 100644
--- a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh
+++ b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh
@@ -449,7 +449,7 @@ private:
//////////////////////////////////////
priVars[energyEq0Idx] = thisTemperature;
priVars[energyEqSolidIdx] = thisTemperature;
- Scalar capPress[numPhases];
+ std::vector<Scalar> capPress(numPhases);
//obtain pc according to saturation
const auto &materialParams =