Commit e5a9b993 authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[freeflow]

updated call of masstransfermodel class

[multidomain]
- if a boundarylayer model for mass transfer is used, now this also means a boundary layer
  model for energy transfer is used too (which is more consistent)
- added a test for the masstransfer and boundary layer models

reviewed by gruenich



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14757 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent d07c495c
......@@ -44,19 +44,19 @@ class MassTransferModel
*
* \param saturation Saturation of the component's phase inside the porous medium
* \param porosity Porosity of the porous medium
* \param charPoreRadius Characteristic pore radius of the porous medium
* \param blThickness Free flow boundary layer thickness
* \param massTransferModel The index of the chosen mass transfer model
*/
public:
MassTransferModel(Scalar saturation, Scalar porosity, Scalar charPoreRadius,
MassTransferModel(Scalar saturation, Scalar porosity,
Scalar blThickness, unsigned int massTransferModel)
: saturation_(saturation), porosity_(porosity), charPoreRadius_(charPoreRadius),
: saturation_(saturation), porosity_(porosity),
blThickness_(blThickness), massTransferModel_(massTransferModel)
{
moistureContent_ = saturation_ * porosity_;
massTransferCoeff_ = 0.0; // dummy default
charPoreRadius_ = 1e10; // dummy default
capillaryPressure_ = 1e10; // dummy default
}
......@@ -64,6 +64,10 @@ public:
void setMassTransferCoeff(Scalar massTransferCoeff)
{ massTransferCoeff_ = massTransferCoeff; }
//! \brief Sets the characteristic pore diameter \f$[m]\f$.
void setCharPoreRadius(Scalar charPoreRadius)
{ charPoreRadius_ = charPoreRadius; }
//! \brief Sets the capillary pressure \f$[Pa]\f$.
void setCapillaryPressure(Scalar capillaryPressure)
{ capillaryPressure_ = capillaryPressure; }
......@@ -98,8 +102,9 @@ public:
// Schlünder model (Schlünder, CES 1988)
else if (massTransferModel_ == 2)
{
Scalar charPoreRadius = charPoreRadius_;
Scalar massTransferCoeff = 1. + 2./M_PI * charPoreRadius / blThickness_
// check if characteristic pore radius was set
assert (charPoreRadius_ < 9e9);
Scalar massTransferCoeff = 1. + 2./M_PI * charPoreRadius_ / blThickness_
* std::sqrt(M_PI/(4.*moistureContent_))
* (std::sqrt(M_PI/(4.*moistureContent_)) - 1.);
......@@ -122,8 +127,9 @@ public:
// modified Schlünder model
else if (massTransferModel_ == 4)
{
Scalar charPoreRadius = charPoreRadius_;
Scalar massTransferCoeff = 1. + 2./M_PI * charPoreRadius / blThickness_
// check if characteristic pore radius was set
assert (charPoreRadius_ < 9e9);
Scalar massTransferCoeff = 1. + 2./M_PI * charPoreRadius_ / blThickness_
* (1./moistureContent_) * (1./moistureContent_ - 1.);
return 1./massTransferCoeff;
......@@ -138,11 +144,11 @@ private:
Scalar saturation_;
Scalar porosity_;
Scalar moistureContent_;
Scalar charPoreRadius_;
Scalar blThickness_;
unsigned int massTransferModel_;
Scalar massTransferCoeff_;
Scalar charPoreRadius_;
Scalar capillaryPressure_;
};
......
......@@ -73,6 +73,11 @@ public:
wPhaseIdx2 = TwoPTwoCNIIndices::wPhaseIdx, //!< Index for the liquid phase
nPhaseIdx2 = TwoPTwoCNIIndices::nPhaseIdx //!< Index for the gas phase
};
enum { nPhaseIdx1 = Stokes2cniIndices::phaseIdx }; //!< Index of the free-flow phase of the fluidsystem
enum { // indices of the components
transportCompIdx1 = Stokes2cniIndices::transportCompIdx, //!< Index of transported component
phaseCompIdx1 = Stokes2cniIndices::phaseCompIdx //!< Index of main component of the phase
};
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, dim> DimVector; //!< A field vector with dim entries
......@@ -112,7 +117,81 @@ public:
if (cParams.boundaryTypes2.isCouplingInflow(energyEqIdx2))
{
if (globalProblem.sdProblem1().isCornerPoint(globalPos1))
unsigned int blModel = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, BoundaryLayer, Model);
// only enter here, if a boundary layer model is used for the computation of the diffusive fluxes
if (blModel)
{
// convective energy flux
const Scalar convectiveFlux =
normalMassFlux1 *
cParams.elemVolVarsCur1[vertInElem1].enthalpy();
// diffusive transported energy
const Scalar massFractionOut = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefMassfrac);
const Scalar M1 = FluidSystem::molarMass(transportCompIdx1);
const Scalar M2 = FluidSystem::molarMass(phaseCompIdx1);
const Scalar X2 = 1.0 - massFractionOut;
const Scalar massToMoleDenominator = M2 + X2*(M1 - M2);
const Scalar moleFractionOut = massFractionOut * M2 /massToMoleDenominator;
Scalar normalMoleFracGrad =
cParams.elemVolVarsCur1[vertInElem1].moleFraction(transportCompIdx1) -
moleFractionOut;
const Scalar velocity = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefVelocity);
// current position + additional virtual runup distance
const Scalar distance = globalPos1[0] + GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, BoundaryLayer, Offset);
const Scalar kinematicViscosity = cParams.elemVolVarsCur1[vertInElem2].kinematicViscosity();
BoundaryLayerModel<TypeTag> boundaryLayerModel(velocity, distance, kinematicViscosity, blModel);
normalMoleFracGrad /= boundaryLayerModel.massBoundaryLayerThickness();
Scalar diffusiveFlux =
bfNormal1.two_norm() *
normalMoleFracGrad *
(boundaryVars1.diffusionCoeff(transportCompIdx1) + boundaryVars1.eddyDiffusivity()) *
boundaryVars1.molarDensity();
// multiply the diffusive flux with the mass transfer coefficient
unsigned int mtModel = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, MassTransfer, Model);
if (mtModel != 0)
{
MassTransferModel<TypeTag> massTransferModel(cParams.elemVolVarsCur2[vertInElem2].saturation(wPhaseIdx2),
cParams.elemVolVarsCur2[vertInElem2].porosity(),
boundaryLayerModel.massBoundaryLayerThickness(),
mtModel);
if (mtModel == 1)
massTransferModel.setMassTransferCoeff(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, MassTransfer, Coefficient));
if (mtModel == 2 || mtModel == 4)
massTransferModel.setCharPoreRadius(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, MassTransfer, CharPoreRadius));
if (mtModel == 3)
massTransferModel.setCapillaryPressure(cParams.elemVolVarsCur2[vertInElem2].capillaryPressure());
diffusiveFlux *= massTransferModel.massTransferCoefficient();
}
// the boundary layer models only work for two components
assert(numComponents == 2);
Scalar diffusiveEnergyFlux = 0.0;
diffusiveEnergyFlux += diffusiveFlux * FluidSystem::molarMass(transportCompIdx1)
* boundaryVars1.componentEnthalpy(transportCompIdx1);
diffusiveEnergyFlux -= diffusiveFlux * FluidSystem::molarMass(phaseCompIdx1)
* boundaryVars1.componentEnthalpy(phaseCompIdx1);
// conductive transported energy
const Scalar temperatureOut = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefTemperature);
Scalar normalTemperatureGrad =
cParams.elemVolVarsCur1[vertInElem1].temperature() -
temperatureOut;
normalTemperatureGrad /= boundaryLayerModel.thermalBoundaryLayerThickness();
const Scalar conductiveFlux =
bfNormal1.two_norm() *
normalTemperatureGrad *
(boundaryVars1.thermalConductivity() + boundaryVars1.thermalEddyConductivity());
couplingRes2.accumulate(lfsu_n.child(energyEqIdx2), vertInElem2,
-(convectiveFlux - diffusiveEnergyFlux - conductiveFlux));
}
else if (globalProblem.sdProblem1().isCornerPoint(globalPos1))
{
const Scalar convectiveFlux =
normalMassFlux1 *
......
......@@ -401,11 +401,12 @@ class TwoCStokesTwoPTwoCLocalOperator :
{
MassTransferModel<TypeTag> massTransferModel(cParams.elemVolVarsCur2[vertInElem2].saturation(wPhaseIdx2),
cParams.elemVolVarsCur2[vertInElem2].porosity(),
GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, MassTransfer, CharPoreDiameter),
boundaryLayerModel.massBoundaryLayerThickness(),
massTransferModel_);
if (massTransferModel_ == 1)
massTransferModel.setMassTransferCoeff(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, MassTransfer, Coefficient));
if (massTransferModel_ == 2 || massTransferModel_ == 4)
massTransferModel.setCharPoreRadius(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, MassTransfer, CharPoreRadius));
if (massTransferModel_ == 3)
massTransferModel.setCapillaryPressure(cParams.elemVolVarsCur2[vertInElem2].capillaryPressure());
const Scalar massTransferCoeff = massTransferModel.massTransferCoefficient();
......
......@@ -17,3 +17,23 @@ add_dumux_test(test_2cnistokes2p2cni_2pm test_2cnistokes2p2cni test_2cnistokes2p
${CMAKE_CURRENT_BINARY_DIR}/test_2cnistokes2p2cni
-ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/test_2cnistokes2p2cni_reference.input
)
# test: test_2cnistokes2p2cni_boundarylayer_1ff
add_dumux_test(test_2cnistokes2p2cni_boundarylayer_1ff test_2cnistokes2p2cni test_2cnistokes2p2cni.cc
${CMAKE_SOURCE_DIR}/bin/runTest.sh
${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
${CMAKE_SOURCE_DIR}/test/references/2cnistokes2p2cniboundarylayer-ff-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/stokes2cni-00004.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_2cnistokes2p2cni
-ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/test_2cnistokes2p2cni_boundarylayer.input
)
# test: test_2cnistokes2p2cni_boundarylayer_2pm
add_dumux_test(test_2cnistokes2p2cni_boundarylayer_2pm test_2cnistokes2p2cni test_2cnistokes2p2cni.cc
${CMAKE_SOURCE_DIR}/bin/runTest.sh
${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
${CMAKE_SOURCE_DIR}/test/references/2cnistokes2p2cniboundarylayer-pm-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/2p2cni-00004.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_2cnistokes2p2cni
-ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/test_2cnistokes2p2cni_boundarylayer.input
)
......@@ -58,7 +58,7 @@ void printUsage(const char *progName, const std::string &errorMsg)
"\n"
"[MassTransferModel]\n"
"Coefficient Coeffient used for the exponential law (model 1) [-]\n"
"CharPoreDiameter Characteristic pore diameter for Schluender model (model 2+4) [m]\n"
"CharPoreRadius Characteristic pore radius for Schluender model (model 2+4) [m]\n"
"\n";
std::cout << errorMessageOut
......
[TimeManager]
DtInitial = 5e-1 # [s]
MaxTimeStepSize = 360 # [s]
InitTime = 0 # [s] Initialization time without coupling
TEnd= 3600 # [s]
EpisodeLength = 43200 # [s]
[Grid]
UseInterfaceMeshCreator = true
# Grid details
File = grids/interfacedomain.dgf
XMin = 0.0 # [m]
XMax = 0.25 # [m]
YMin = 0.0 # [m]
YMax = 0.5 # [m]
RunUpDistanceX = 0.0 # [m] Horizontal position without coupling to PM
NoDarcyX = 0.0 # [m] Horizontal position without PM below
InterfacePosY = 0.25 # [m] Vertical position of coupling interface
# Number of elements in x-, y-, z-direction
CellsX = 15
CellsY = 20
# Grading and refinement of the mesh in y direction
GradingFactorY = 1.1
Refinement = 0
[Output]
NameFF = stokes2cni
NamePM = 2p2cni
#Frequency of restart file, flux and VTK output
FreqRestart = 1000 # how often restart files are written out
FreqOutput = 10 # frequency of VTK output
FreqMassOutput = 2 # frequency of mass and evaporation rate output (Darcy)
FreqFluxOutput = 1000 # frequency of detailed flux output
FreqVaporFluxOutput = 2 # frequency of summarized flux output
LiveEvaporationRates = false # plot evaporation rates using gnuplot interface
[Stokes]
StabilizationAlpha = -1.0
[FreeFlow]
UseDirichletBC = true # dirichlet values are set at the inflow boundary
RefVelocity = 3.5 # [m/s]
RefPressure = 1e5 # [Pa]
RefMassfrac = 0.008 # [-]
RefTemperature = 298.15 # [K]
SinusVelAmplitude = 0.0 # [m/s]
SinusVelPeriod = 3600 # [s]
SinusPressureAmplitude = 0.0 # [Pa]
SinusPressurePeriod = 3600 # [s]
SinusConcentrationAmplitude = 0.0 # [-]
SinusConcentrationPeriod = 3600 # [s]
SinusTemperatureAmplitude = 0.0 # [K]
SinusTemperaturePeriod = 3600 # [s]
[BoundaryLayer]
Model = 4 # Number/ID of the used model
Offset = 0.25 # Virtual run-up distance for BL models [m]
YPlus = 10 # Conversion value (model 4-6) [-]
[MassTransfer]
Model = 4
CharPoreRadius = 0.001 # Characteristic pore radius for Schluender model (model 2+4) [m]
[PorousMedium]
RefPressure = 1e5 # [Pa]
RefTemperature = 298.15 # [K]
InitialSw = 0.98 # [-]
[SpatialParams]
AlphaBJ = 1.0 # [-]
Permeability = 2.65e-10 # [m^2]
Porosity = 0.41 # [-]
Swr = 0.005 # [-]
Snr = 0.01 # [-]
VgAlpha = 6.371e-4 # [1/Pa]
VgN = 8.0 # [-]
Pentry = 1357
BCLambda = 6.960
LambdaSolid = 5.26 # [W/(m*K)]
[Newton]
MaxRelativeShift = 1e-5
TargetSteps = 8
MaxSteps = 12
WriteConvergence = false
MaxTimeStepDivisions = 20
[LinearSolver]
ResidualReduction = 1e-9
Verbosity = 0
MaxIterations = 200
\ No newline at end of file
......@@ -58,7 +58,7 @@ void printUsage(const char *progName, const std::string &errorMsg)
"\n"
"[MassTransferModel]\n"
"Coefficient Coeffient used for the exponential law (model 1) [-]\n"
"CharPoreDiameter Characteristic pore diameter for Schluender model (model 2+4) [m]\n"
"CharPoreRadius Characteristic pore radius for Schluender model (model 2+4) [m]\n"
"\n";
std::cout << errorMessageOut
......
This diff is collapsed.
This diff is collapsed.
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