Commit e052e433 authored by Philipp Nuske's avatar Philipp Nuske
Browse files

Choosing which parametrization is used for the Nusselt / Sherwood correlation:

Introducing a new property:
NusseltFormulation
SherwoodFormulation

Defining possible values in dimensionlessnumbers

The actual call to the function now additionally needs a flag for which formulation to choose. 

add a bunch of consts in order not to accidentially change parameteres (mpncvoluemvariablesiakinetic)

reviewed by Bernd

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@11504 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 8030cc7e
......@@ -95,6 +95,17 @@ static Scalar prandtlNumber(const Scalar dynamicViscosity,
return dynamicViscosity * heatCapacity / thermalConductivity;
}
/*!
* \brief A container for possible values of the property for selecting which nusselt parametrization to choose.
* The actual value is set vie the property NusseltFormulation
*/
struct NusseltFormulation
{
static const int dittusBoelter = 0;
static const int WakaoKaguei = 1;
static const int VDI = 2;
};
/*!
* \brief Calculate the Nusselt Number [-] (Nu).
*
......@@ -116,41 +127,51 @@ static Scalar prandtlNumber(const Scalar dynamicViscosity,
* \param reynoldsNumber Dimensionless number relating inertial and viscous forces [-].
* \param prandtlNumber Dimensionless number relating viscosity and thermal diffusivity (temperaturleitfaehigkeit) [-].
* \param porosity The fraction of the porous medium which is void space.
*
* \param formulation Switch for deciding which parametrization of the Nusselt number is to be used. Set via the property NusseltFormulation.
* \return The Nusselt number as calculated from the input parameters [-].
*/
static Scalar nusseltNumberForced(const Scalar reynoldsNumber,
const Scalar prandtlNumber,
const Scalar porosity)
const Scalar porosity,
const int formulation )
{
/* example: very common and simple case: flow straight circular pipe, only convection (no boiling), 10000<Re<120000, 0.7<Pr<120, far from pipe entrance, smooth surface of pipe ...
Dittus, F.W and Boelter, L.M.K, Heat Transfer in Automobile Radiators of the Tubular Type, Publications in Engineering, Vol. 2, pages 443-461, 1930
*/
// return 0.023 * pow(reynoldsNumber, 0.8) * pow(prandtlNumber,0.33);
/* example: flow through porous medium *single phase*, fit to many different data
* Wakao and Kaguei, Heat and mass Transfer in Packed Beds, Gordon and Breach Science Publishers, page 293
*/
return 2. + 1.1 * pow(prandtlNumber,(1./3.)) * pow(reynoldsNumber, 0.6);
/* example: VDI Waermeatlas 10. Auflage 2006, flow in packed beds, page Gj1, see also other sources and limitations therein.
* valid for 0.1<Re<10000, 0.6<Pr/Sc<10000, packed beds of perfect spheres.
*
*/
// Scalar numerator = 0.037 * pow(reynoldsNumber,0.8) * prandtlNumber ;
// Scalar reToMin01 = pow(reynoldsNumber,-0.1);
// Scalar prTo23 = pow(prandtlNumber, (2./3. ) ) ; // MIND THE pts! :-( otherwise the integer exponent version is chosen
// Scalar denominator = 1+ 2.443 * reToMin01 * (prTo23 -1.) ;
//
// Scalar nusseltTurbular = numerator / denominator;
// Scalar nusseltLaminar = 0.664 * sqrt(reynoldsNumber) * pow(prandtlNumber, (1./3.) );
// Scalar nusseltSingleSphere = 2 + sqrt( pow(nusseltLaminar,2.) + pow(nusseltTurbular,2.));
//
// Scalar funckyFactor = 1 + 1.5 * (1.-porosity); // for spheres of same size
// Scalar nusseltNumber = funckyFactor * nusseltSingleSphere ;
//
// return nusseltNumber;
if (formulation == NusseltFormulation::dittusBoelter){
/* example: very common and simple case: flow straight circular pipe, only convection (no boiling), 10000<Re<120000, 0.7<Pr<120, far from pipe entrance, smooth surface of pipe ...
Dittus, F.W and Boelter, L.M.K, Heat Transfer in Automobile Radiators of the Tubular Type, Publications in Engineering, Vol. 2, pages 443-461, 1930
*/
return 0.023 * pow(reynoldsNumber, 0.8) * pow(prandtlNumber,0.33);
}
else if (formulation == NusseltFormulation::WakaoKaguei){
/* example: flow through porous medium *single phase*, fit to many different data
* Wakao and Kaguei, Heat and mass Transfer in Packed Beds, Gordon and Breach Science Publishers, page 293
*/
return 2. + 1.1 * pow(prandtlNumber,(1./3.)) * pow(reynoldsNumber, 0.6);
}
else if (formulation == NusseltFormulation::VDI){
/* example: VDI Waermeatlas 10. Auflage 2006, flow in packed beds, page Gj1, see also other sources and limitations therein.
* valid for 0.1<Re<10000, 0.6<Pr/Sc<10000, packed beds of perfect spheres.
*
*/
Scalar numerator = 0.037 * pow(reynoldsNumber,0.8) * prandtlNumber ;
Scalar reToMin01 = pow(reynoldsNumber,-0.1);
Scalar prTo23 = pow(prandtlNumber, (2./3. ) ) ; // MIND THE pts! :-( otherwise the integer exponent version is chosen
Scalar denominator = 1+ 2.443 * reToMin01 * (prTo23 -1.) ;
Scalar nusseltTurbular = numerator / denominator;
Scalar nusseltLaminar = 0.664 * sqrt(reynoldsNumber) * pow(prandtlNumber, (1./3.) );
Scalar nusseltSingleSphere = 2 + sqrt( pow(nusseltLaminar,2.) + pow(nusseltTurbular,2.));
Scalar funckyFactor = 1 + 1.5 * (1.-porosity); // for spheres of same size
Scalar nusseltNumber = funckyFactor * nusseltSingleSphere ;
return nusseltNumber;
}
else {
DUNE_THROW(Dune::NotImplemented, "wrong index");
}
}
......@@ -180,6 +201,15 @@ static Scalar schmidtNumber(const Scalar dynamicViscosity,
return dynamicViscosity / (massDensity * diffusionCoefficient);
}
/*!
* \brief A container for possible values of the property for selecting which sherwood parametrization to choose.
* The actual value is set vie the property SherwoodFormulation
*/
struct SherwoodFormulation
{
static const int WakaoKaguei = 0;
};
/*!
* \brief Calculate the Sherwood Number [-] (Sh).
*
......@@ -204,17 +234,24 @@ static Scalar schmidtNumber(const Scalar dynamicViscosity,
*
* \param schmidtNumber Dimensionless number relating viscosity and mass diffusivity [-].
* \param reynoldsNumber Dimensionless number relating inertial and viscous forces [-].
* \param formulation Switch for deciding which parametrization of the Sherwood number is to be used. Set via the property SherwoodFormulation.
* \return The Nusselt number as calculated from the input parameters [-].
*/
static Scalar sherwoodNumber(const Scalar reynoldsNumber,
const Scalar schmidtNumber)
const Scalar schmidtNumber,
const int formulation)
{
/* example: flow through porous medium *single phase*
* Wakao and Kaguei, Heat and mass Transfer in Packed Beds, Gordon and Breach Science Publishers, page 156
*/
return 2. + 1.1 * pow(schmidtNumber,(1./3.)) * pow(reynoldsNumber, 0.6);
}
if (formulation == SherwoodFormulation::WakaoKaguei){
/* example: flow through porous medium *single phase*
* Wakao and Kaguei, Heat and mass Transfer in Packed Beds, Gordon and Breach Science Publishers, page 156
*/
return 2. + 1.1 * pow(schmidtNumber,(1./3.)) * pow(reynoldsNumber, 0.6);
}
else {
DUNE_THROW(Dune::NotImplemented, "wrong index");
}}
/*!
......
......@@ -106,7 +106,9 @@ public:
if (phaseIdx not_eq sPhaseIdx) {
storage[energyEq0Idx + phaseIdx] +=
fs.density(phaseIdx) *
fs.internalEnergy(phaseIdx) *
// fs.internalEnergy(phaseIdx) *
#warning balancing enthalpy (no volume change work)
fs.enthalpy(phaseIdx) *
fs.saturation(phaseIdx) *
volVars.porosity();
}
......
......@@ -116,34 +116,6 @@ public:
void checkDefined() const
{
}
/*!
* \brief Check the set variables as to whether they are in physically possible ranges.
*
* \param fluidState Container for all the secondary variables concerning the fluids
* \param globalPos The position at which the check is conducted
*
* Since we are isothermal, we don't need to do anything!
*/
bool physicalness(const FluidState & fluidState,
const GlobalPosition & globalPos)
{
return true; // all the checks went through: tell calling function, nothing bad could be found.
}
/*!
* \brief Output for the case that the current state is not physical.
* This calls the output functions of the modules and throws and exception:
* i.e. a smaller timestep is tried.
*
* Since we are isothermal, we don't need to do anything!
*
* \param fs Container for all the secondary variables concerning the fluids
* \param message A string returning the error message for this module
*/
const void physicalnessError(const FluidState & fs,
std::stringstream & message)
{ }
};
/*!
......@@ -257,38 +229,6 @@ public:
Valgrind::CheckDefined(soilDensity_);
}
/*!
* \brief Check whether the calculated values are reasonable.
*
* \param fs Container for all the secondary variables concerning the fluids
* \param globalPos The position at which the check is conducted
*/
bool physicalness(const FluidState & fs,
const GlobalPosition & globalPos)
{
const Scalar eps = 1e-6 ;
const Scalar temperatureTest = fs.temperature(/*dummy=*/0);
if (not std::isfinite(temperatureTest)
or temperatureTest < 0.-eps )
return false; // unphysical value found: tell calling function, sth went wrong!
return true; // all the checks went through: tell calling function, nothing bad could be found.
}
/*!
* \brief Output for the case that the current state is not physical.
* This is called if the physicalness funcitons returned false.
*
* \param fs Container for all the secondary variables concerning the fluids
* \param message A string returning the error message for this module
*/
const void physicalnessError(const FluidState & fs,
std::stringstream & message)
{
message <<"Energy: \n";
for(int energyEqIdx=0; energyEqIdx<numEnergyEqs; energyEqIdx++)
message << "\tT" <<"_"<<FluidSystem::phaseName(energyEqIdx)<<"="<< fs.temperature(energyEqIdx) <<"\n";
}
protected:
Scalar heatCapacity_;
Scalar soilDensity_;
......
......@@ -116,46 +116,6 @@ public:
}
}
/*!
* \brief Check whether the calculated values are reasonable.
*
* \param fs Container for all the secondary variables concerning the fluids
* \param globalPos The position at which the check is conducted
*/
bool physicalness(const FluidState & fs,
const GlobalPosition & globalPos)
{
const Scalar eps = 1e-6 ;
for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++){
for (int compIdx=0; compIdx< numComponents; ++ compIdx){
const Scalar xTest = fs.moleFraction(phaseIdx, compIdx);
if (not std::isfinite(xTest)
or xTest < 0.-eps
or xTest > 1.+eps )
return false; // unphysical value found: tell calling function, sth went wrong!
}
}
return true; // all the checks went through: tell calling function, nothing bad could be found.
}
/*!
* \brief Output for the case that the current state is not physical.
* This is called if the physicalness funcitons returned false.
*
* \param fs Container for all the secondary variables concerning the fluids
* \param message A string returning the error message for this module
*/
const void physicalnessError(const FluidState & fs,
std::stringstream & message)
{
message <<"Mass: \n";
for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++)
for (int compIdx=0; compIdx< numComponents; ++ compIdx)
message << "\tx" <<"_"<<FluidSystem::phaseName(phaseIdx)
<<"^"<<FluidSystem::componentName(compIdx)<<"="
<< fs.moleFraction(phaseIdx, compIdx) <<"\n";
}
/*!
* \brief If running in valgrind this makes sure that all
* quantities in the volume variables are defined.
......
......@@ -62,6 +62,12 @@ NEW_PROP_TAG(EnableKineticEnergy);
//! average the velocity in the model
NEW_PROP_TAG(VelocityAveragingInModel);
//! which nusselt correlation to use
NEW_PROP_TAG(NusseltFormulation);
//! which sherwood correlation to use
NEW_PROP_TAG(SherwoodFormulation);
}
}
......
......@@ -45,7 +45,6 @@
#include "energy/mpncvolumevariablesenergykinetic.hh"
#include "energy/mpncvtkwriterenergykinetic.hh"
namespace Dumux
{
namespace Properties
......@@ -56,48 +55,9 @@ namespace Properties
*/
SET_TYPE_PROP(BoxMPNCKinetic, Model, MPNCModelKinetic<TypeTag>);
/*!
* \brief Set the property for the awn surface by retrieving it from
* the spatial parameters.
*/
SET_PROP(BoxMPNCKinetic, AwnSurface)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
public:
typedef typename SpatialParams::AwnSurface type;
};
/*!
* \brief Set the property for the awn surface by retrieving it from
* the spatial parameters.
*/
SET_PROP(BoxMPNCKinetic, AnsSurface)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
public:
typedef typename SpatialParams::AnsSurface type;
};
/*!
* \brief Set the property for the awn surface by retrieving it from
* the spatial parameters.
*/
SET_PROP(BoxMPNCKinetic, AwsSurface)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
public:
typedef typename SpatialParams::AwsSurface type;
};
/*!
* \brief Set the property for the awn surface parameters by extracting
* it from awn surface??.
* it from awn surface.
*/
SET_PROP(BoxMPNCKinetic, AwnSurfaceParams)
{
......@@ -109,8 +69,8 @@ public:
};
/*!
* \brief Set the property for the awn surface parameters by extracting
* it from awn surface??.
* \brief Set the property for the aws surface parameters by extracting
* it from aws surface.
*/
SET_PROP(BoxMPNCKinetic, AwsSurfaceParams)
{
......@@ -122,8 +82,8 @@ public:
};
/*!
* \brief Set the property for the awn surface parameters by extracting
* it from awn surface??.
* \brief Set the property for the ans surface parameters by extracting
* it from the surface.
*/
SET_PROP(BoxMPNCKinetic, AnsSurfaceParams)
{
......@@ -135,7 +95,30 @@ public:
};
SET_BOOL_PROP(BoxMPNCKinetic, VelocityAveragingInModel, true);
}
}
/*!
* \brief Set the default formulation for the nusselt correlation
* Other possible parametrizations can be found in the dimensionlessnumbers
*/
SET_PROP(BoxMPNCKinetic, NusseltFormulation ){
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef DimensionlessNumbers<Scalar> DimLessNum;
public:
static constexpr int value = DimLessNum::NusseltFormulation::WakaoKaguei;};
/*!
* \brief Set the default formulation for the sherwood correlation
* Other possible parametrizations can be found in the dimensionlessnumbers
*/
SET_PROP(BoxMPNCKinetic, SherwoodFormulation ){
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef DimensionlessNumbers<Scalar> DimLessNum;
public:
static constexpr int value = DimLessNum::SherwoodFormulation::WakaoKaguei;};
} // end namespace Properties
} // end namespace Dumux
#endif
......@@ -54,7 +54,10 @@ class MPNCVolumeVariablesIA<TypeTag, enableKinetic, /*bool enableKineticEnergy=*
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
......@@ -64,8 +67,12 @@ class MPNCVolumeVariablesIA<TypeTag, enableKinetic, /*bool enableKineticEnergy=*
enum { nCompIdx = FluidSystem::nCompIdx } ;
enum { wCompIdx = FluidSystem::wCompIdx } ;
enum { dim = GridView::dimension};
enum { numEnergyEqs = Indices::NumPrimaryEnergyVars};
enum { nusseltFormulation = GET_PROP_VALUE(TypeTag, NusseltFormulation)} ;
enum { sherwoodFormulation = GET_PROP_VALUE(TypeTag, SherwoodFormulation)} ;
typedef DimensionlessNumbers<Scalar> DimLessNum ;
typedef DimensionlessNumbers<Scalar> DimLessNum;
typedef Dune::FieldVector<Scalar,dim> GlobalPosition;
......@@ -82,16 +89,6 @@ class MPNCVolumeVariablesIA<TypeTag, enableKinetic, /*bool enableKineticEnergy=*
public:
/*!
* \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
*
* \param volVars The volume variables
* \param fluisState Container for all the secondary variables concerning the fluids
* \param paramCache Container for cache parameters
* \param priVars The primary Variables
* \param problem The problem
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the fully implicit scheme
* \param scvIdx The index of the sub-control volumete element
*
*/
void update(const VolumeVariables & volVars,
const FluidState & fluidState,
......@@ -102,14 +99,17 @@ public:
const FVElementGeometry & fvGeometry,
const unsigned int scvIdx)
{
// obtain (standard) material parameters (needed for the residual saturations)
const MaterialLawParams & materialParams = problem.spatialParams().materialLawParams(element,fvGeometry,scvIdx) ;
//obtain parameters for interfacial area constitutive relations
AwnSurfaceParams aWettingNonWettingSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element,fvGeometry,scvIdx) ;
AnsSurfaceParams aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element,fvGeometry,scvIdx) ;
const AwnSurfaceParams & aWettingNonWettingSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element,fvGeometry,scvIdx) ;
const AnsSurfaceParams & aNonWettingSolidSurfaceParams = problem.spatialParams().aNonWettingSolidSurfaceParams(element,fvGeometry,scvIdx) ;
Valgrind::CheckDefined(aWettingNonWettingSurfaceParams);
Valgrind::CheckDefined(aNonWettingSolidSurfaceParams);
const Scalar pc = fluidState.pressure(nPhaseIdx) - fluidState.pressure(wPhaseIdx);
const Scalar pc = fluidState.pressure(nPhaseIdx) - fluidState.pressure(wPhaseIdx) ;
const Scalar Sw = fluidState.saturation(wPhaseIdx) ;
Valgrind::CheckDefined(Sw);
Valgrind::CheckDefined(pc);
......@@ -134,25 +134,31 @@ public:
else
#endif
awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, Sw, pc ); // 10.; //
if (AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc ) < 0)
{
Scalar dummy = 0 ;
}
awn = AwnSurface::interfacialArea(aWettingNonWettingSurfaceParams, materialParams, Sw, pc ); // 10.; //
interfacialArea_[wPhaseIdx][nPhaseIdx] = awn ; //10. ;//
interfacialArea_[nPhaseIdx][wPhaseIdx] = interfacialArea_[wPhaseIdx][nPhaseIdx];
interfacialArea_[wPhaseIdx][wPhaseIdx] = 0. ;
Scalar ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, Sw, pc ); // 10.; //
Scalar ans = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams,Sw, pc ); // 10.; //
// if (ans <0 )
// ans = 0 ;
// Switch for using a a_{wn} relations that has some "maximum capillary pressure" as parameter.
// That value is obtained by regularization of the pc(Sw) function.
#if USE_PCMAX
const Scalar pcMax = problem.spatialParams().pcMax(element,
const Scalar pcMax = problem.spatialParams().pcMax(element,
fvGeometry,
scvIdx);
// I know the solid surface from the pore network. But it is more consistent to use the fit value.
// solidSurface_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.soil.specificSolidsurface);
solidSurface_ = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, /*Sw=*/0., pcMax );
solidSurface_ = AnsSurface::interfacialArea(aNonWettingSolidSurfaceParams, materialParams, /*Sw=*/0., pcMax );
Valgrind::CheckDefined(solidSurface_);
// if (ans > solidSurface_){
......@@ -181,9 +187,9 @@ public:
interfacialArea_[sPhaseIdx][wPhaseIdx] = interfacialArea_[wPhaseIdx][sPhaseIdx];
interfacialArea_[sPhaseIdx][sPhaseIdx] = 0. ;
#else
AwsSurfaceParams aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams();
const AwsSurfaceParams & aWettingSolidSurfaceParams = problem.spatialParams().aWettingSolidSurfaceParams();
Valgrind::CheckDefined(aWettingSolidSurfaceParams);
const Scalar aws = AwsSurface::interfacialArea(aWettingSolidSurfaceParams, Sw, pc ); // 10.; //
const Scalar aws = AwsSurface::interfacialArea(aWettingSolidSurfaceParams,materialParams, Sw, pc ); // 10.; //
interfacialArea_[wPhaseIdx][sPhaseIdx] = aws ;
interfacialArea_[sPhaseIdx][wPhaseIdx] = interfacialArea_[wPhaseIdx][sPhaseIdx];
interfacialArea_[sPhaseIdx][sPhaseIdx] = 0. ;
......@@ -235,14 +241,16 @@ public:
nusseltNumber_[phaseIdx] = DimLessNum::nusseltNumberForced(reynoldsNumber_[phaseIdx],
prandtlNumber_[phaseIdx],
porosity);
porosity,
nusseltFormulation);
schmidtNumber_[phaseIdx] = DimLessNum::schmidtNumber(dynamicViscosity,
density,
diffCoeff);
sherwoodNumber_[phaseIdx] = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
schmidtNumber_[phaseIdx]);
schmidtNumber_[phaseIdx],
sherwoodFormulation);
}
}
......@@ -251,9 +259,6 @@ public:
*
* This is _only_ required by the kinetic mass/energy modules
*
* \param phaseIIdx The local index of the first phase
* \param phaseJIdx The local index of the second phase
*
*/
const Scalar interfacialArea(const unsigned int phaseIIdx, const unsigned int phaseJIdx) const
{
......@@ -356,20 +361,14 @@ class MPNCVolumeVariablesIA<TypeTag, /*enableKinetic=*/true, /*bool enableKineti
enum { nCompIdx = FluidSystem::nCompIdx };
enum { dim = GridView::dimension};
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
enum { sherwoodFormulation = GET_PROP_VALUE(TypeTag, SherwoodFormulation)} ;
typedef Dune::FieldVector<Scalar,dim> GlobalPosition;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
public:
/*!
* \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
*
* \param volVars The volume variables
* \param fluisState Container for all the secondary variables concerning the fluids
* \param paramCache Container for cache parameters
* \param priVars The primary Variables
* \param problem The problem
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the fully implicit scheme
* \param scvIdx The index of the sub-control volumete element
*/
void update(const VolumeVariables & volVars,
const FluidState & fluidState,
......@@ -381,14 +380,17 @@ public:
const unsigned int scvIdx)
{
//obtain parameters for awnsurface
AwnSurfaceParams awnSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element,fvGeometry,scvIdx) ;
const AwnSurfaceParams & awnSurfaceParams = problem.spatialParams().aWettingNonWettingSurfaceParams(element,fvGeometry,scvIdx) ;
// obtain (standard) material parameters (needed for the residual saturations)
const MaterialLawParams & materialParams = problem.spatialParams().materialLawParams(element,fvGeometry,scvIdx) ;
Valgrind::CheckDefined(awnSurfaceParams);
const Scalar Sw = fluidState.saturation(wPhaseIdx) ;
const Scalar pc = fluidState.pressure(nPhaseIdx) - fluidState.pressure(wPhaseIdx);
// so far there is only a model for kinetic mass transfer between fluid phases
interfacialArea_ = AwnSurface::interfacialArea(awnSurfaceParams, Sw, pc );
interfacialArea_ = AwnSurface::interfacialArea(awnSurfaceParams, materialParams, Sw, pc );
Valgrind::CheckDefined(interfacialArea_);
......@@ -420,15 +422,13 @@ public:
diffCoeff);
sherwoodNumber_[phaseIdx] = DimLessNum::sherwoodNumber(reynoldsNumber_[phaseIdx],
schmidtNumber_[phaseIdx]);
schmidtNumber_[phaseIdx],
sherwoodFormulation);
}
}
/*!
* \brief The specific interfacial area between two fluid phases [m^2 / m^3]
*
* \param phaseIIdx The local index of the first phase
* \param phaseJIdx The local index of the second phase
*/
const Scalar interfacialArea(const unsigned int