35#ifndef DUMUX_EVAPORATION_ATMOSPHERE_PROBLEM_HH
36#define DUMUX_EVAPORATION_ATMOSPHERE_PROBLEM_HH
38#include <dune/grid/yaspgrid.hh>
65template <
class TypeTag>
66class EvaporationAtmosphereProblem;
76template<
class TypeTag>
77struct Grid<TypeTag, TTag::EvaporationAtmosphere> {
using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
80template<
class TypeTag>
84template<
class TypeTag>
92template<
class TypeTag>
100template<
class TypeTag>
101struct Scalar<TypeTag, TTag::EvaporationAtmosphere> {
using type = double; };
104template<
class TypeTag>
113template<
class TypeTag>
129template <
class TypeTag>
141 using SubControlVolume =
typename FVElementGeometry::SubControlVolume;
142 using SubControlVolumeFace =
typename FVElementGeometry::SubControlVolumeFace;
144 using Element =
typename GridView::template Codim<0>::Entity;
145 using GlobalPosition =
typename Element::Geometry::GlobalCoordinate;
149 using ParameterCache =
typename FluidSystem::ParameterCache;
155 enum { dimWorld = GridView::dimensionworld };
156 enum { numPhases = ModelTraits::numFluidPhases() };
157 enum { numComponents = ModelTraits::numFluidComponents() };
158 enum { s0Idx = Indices::s0Idx };
159 enum { p0Idx = Indices::p0Idx };
160 enum { conti00EqIdx = Indices::conti0EqIdx };
161 enum { energyEq0Idx = Indices::energyEqIdx };
162 enum { liquidPhaseIdx = FluidSystem::liquidPhaseIdx };
163 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
164 enum { wCompIdx = FluidSystem::H2OIdx };
165 enum { nCompIdx = FluidSystem::N2Idx };
166 enum { numEnergyEqFluid = ModelTraits::numEnergyEqFluid() };
167 enum { numEnergyEqSolid = ModelTraits::numEnergyEqSolid() };
169 static constexpr bool enableChemicalNonEquilibrium = getPropValue<TypeTag, Properties::EnableChemicalNonEquilibrium>();
173 static constexpr auto pressureFormulation = ModelTraits::pressureFormulation();
181 percentOfEquil_ = getParam<Scalar>(
"BoundaryConditions.percentOfEquil");
182 nTemperature_ = getParam<Scalar>(
"FluidSystem.nTemperature");
183 nPressure_ = getParam<Scalar>(
"FluidSystem.nPressure");
184 outputName_ = getParam<std::string>(
"Problem.Name");
185 TInitial_ = getParam<Scalar>(
"InitialConditions.TInitial");
186 SwPMInitial_ = getParam<Scalar>(
"InitialConditions.SwPMInitial");
187 SwFFInitial_ = getParam<Scalar>(
"InitialConditions.SwFFInitial");
188 pnInitial_ = getParam<Scalar>(
"InitialConditions.pnInitial");
189 pnInjection_ = getParam<Scalar>(
"InitialConditions.pnInjection");
190 TInject_ = getParam<Scalar>(
"BoundaryConditions.TInject");
191 massFluxInjectedPhase_ = getParam<Scalar>(
"BoundaryConditions.massFluxInjectedPhase");
194 FluidSystem::init(TInitial_ - 15.0, 453.15, nTemperature_,
195 0.75*pnInitial_, 2.25*pnInitial_, nPressure_);
202 {
return *gridVariables_; }
215 {
return outputName_ ; }
232 BoundaryTypes bcTypes;
234 bcTypes.setAllNeumann();
237 if(onRightBoundary_(globalPos) && globalPos[1] > this->
spatialParams().heightPM() + eps_)
239 bcTypes.setAllDirichlet();
252 return initial_(globalPos);
259 const FVElementGeometry& fvGeometry,
260 const ElementVolumeVariables& elemVolVars,
261 const ElementFluxVariablesCache& elemFluxVarsCache,
262 const SubControlVolumeFace& scvf)
const
264 NumEqVector values(0.0);
265 const auto& globalPos = fvGeometry.scv(scvf.insideScvIdx()).dofPosition();
266 const Scalar massFluxInjectedPhase = massFluxInjectedPhase_ ;
268 ParameterCache dummyCache;
269 FluidState fluidState;
271 for (
int phaseIdx=0; phaseIdx<numPhases; phaseIdx++)
273 fluidState.setPressure(phaseIdx, pnInitial_);
276 fluidState.setTemperature(gasPhaseIdx, TInject_ );
277 fluidState.setTemperature(liquidPhaseIdx, TInitial_ );
284 fluidState.setMoleFraction(gasPhaseIdx, wCompIdx, fluidState.moleFraction(gasPhaseIdx, wCompIdx)*percentOfEquil_ ) ;
285 fluidState.setMoleFraction(gasPhaseIdx, nCompIdx, 1.-fluidState.moleFraction(gasPhaseIdx, wCompIdx) ) ;
291 fluidState.setDensity(gasPhaseIdx,
density);
297 for(
int phaseIdx=0; phaseIdx<numPhases; phaseIdx++)
299 const Scalar h = FluidSystem::enthalpy(fluidState,
302 fluidState.setEnthalpy(phaseIdx, h);
305 const Scalar molarFlux = massFluxInjectedPhase / fluidState.averageMolarMass(gasPhaseIdx);
308 if (onLeftBoundary_(globalPos) && this->
spatialParams().inFF_(globalPos))
310 values[conti00EqIdx + gasPhaseIdx * numComponents + wCompIdx]
311 = -molarFlux * fluidState.moleFraction(gasPhaseIdx, wCompIdx);
312 values[conti00EqIdx + gasPhaseIdx * numComponents + nCompIdx]
313 = -molarFlux * fluidState.moleFraction(gasPhaseIdx, nCompIdx);
315 values[energyEq0Idx + gasPhaseIdx] = - massFluxInjectedPhase
316 * fluidState.enthalpy(gasPhaseIdx) ;
333 return initial_(globalPos);
347 NumEqVector
source(
const Element &element,
348 const FVElementGeometry& fvGeometry,
349 const ElementVolumeVariables& elemVolVars,
350 const SubControlVolume &scv)
const
352 NumEqVector values(0.0);
359 PrimaryVariables initial_(
const GlobalPosition &globalPos)
const
361 PrimaryVariables priVars(0.0);
362 const Scalar T = TInitial_;
366 S[liquidPhaseIdx] = SwPMInitial_;
367 S[gasPhaseIdx] = 1. - S[liquidPhaseIdx] ;
370 S[liquidPhaseIdx] = SwFFInitial_;
371 S[gasPhaseIdx] = 1. - S[liquidPhaseIdx] ;
374 DUNE_THROW(Dune::InvalidStateException,
375 "You should not be here: x=" << globalPos[0] <<
" y= "<< globalPos[dimWorld-1]);
377 for (
int i = 0; i < numPhases - 1; ++i)
378 priVars[s0Idx + i] = S[i];
381 FluidState equilibriumFluidState;
384 for (
int phaseIdx = 0; phaseIdx < numPhases ; ++phaseIdx)
386 equilibriumFluidState.setSaturation(phaseIdx, S[phaseIdx]);
387 equilibriumFluidState.setTemperature(phaseIdx, TInitial_ );
390 const auto &materialParams =
392 std::vector<Scalar> capPress(numPhases);
394 using MaterialLaw =
typename ParentType::SpatialParams::MaterialLaw;
395 using MPAdapter = MPAdapter<MaterialLaw, numPhases>;
397 const int wPhaseIdx = this->
spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
398 MPAdapter::capillaryPressures(capPress, materialParams, equilibriumFluidState, wPhaseIdx);
404 p[liquidPhaseIdx] = pnInitial_ - abs(capPress[liquidPhaseIdx]);
405 p[gasPhaseIdx] = p[liquidPhaseIdx] + abs(capPress[liquidPhaseIdx]);
408 p[gasPhaseIdx] = pnInitial_ ;
409 p[liquidPhaseIdx] = pnInitial_ ;
412 DUNE_THROW(Dune::InvalidStateException,
"You should not be here: x=" << globalPos[0] <<
" y= "<< globalPos[dimWorld-1]);
414 if(pressureFormulation == mostWettingFirst){
417 priVars[p0Idx] = p[liquidPhaseIdx];
419 else if(pressureFormulation == leastWettingFirst){
422 priVars[p0Idx] = p[gasPhaseIdx];
424 else DUNE_THROW(Dune::InvalidStateException,
"EvaporationAtmosphereProblem does not support the chosen pressure formulation.");
426 for (
int energyEqIdx=0; energyEqIdx< numEnergyEqFluid+numEnergyEqSolid; ++energyEqIdx)
427 priVars[energyEq0Idx + energyEqIdx] = T;
429 for (
int phaseIdx=0; phaseIdx<numPhases; phaseIdx++)
430 equilibriumFluidState.setPressure(phaseIdx, p[phaseIdx]);
433 ParameterCache dummyCache;
437 FluidState dryFluidState(equilibriumFluidState);
439 dryFluidState.setMoleFraction(gasPhaseIdx, wCompIdx, dryFluidState.moleFraction(gasPhaseIdx, wCompIdx) * percentOfEquil_ ) ;
440 dryFluidState.setMoleFraction(gasPhaseIdx, nCompIdx, 1.0-dryFluidState.moleFraction(gasPhaseIdx, wCompIdx) ) ;
445 if(enableChemicalNonEquilibrium){
446 for(
int phaseIdx=0; phaseIdx < numPhases; ++ phaseIdx)
448 for(
int compIdx=0; compIdx <numComponents; ++compIdx){
449 int offset = compIdx + phaseIdx * numComponents ;
452 priVars[conti00EqIdx + offset] = equilibriumFluidState.moleFraction(phaseIdx,compIdx) ;
455 priVars[conti00EqIdx + offset] = dryFluidState.moleFraction(phaseIdx,compIdx) ;
458 DUNE_THROW(Dune::InvalidStateException,
"You should not be here: x=" << globalPos[0] <<
" y= "<< globalPos[dimWorld-1]);
465 const Scalar fugH2O = FluidSystem::H2O::vaporPressure(T) ;
466 const Scalar fugN2 = p[gasPhaseIdx] - fugH2O ;
468 priVars[conti00EqIdx + FluidSystem::N2Idx] = fugN2 ;
469 priVars[conti00EqIdx + FluidSystem::H2OIdx] = fugH2O ;
471 Scalar xl[numComponents];
472 Scalar beta[numComponents];
475 const Scalar satVapPressure = FluidSystem::H2O::vaporPressure(TInitial_);
476 xl[FluidSystem::H2OIdx] = x_[liquidPhaseIdx][wCompIdx];
477 xl[FluidSystem::N2Idx] = x_[liquidPhaseIdx][nCompIdx];
478 beta[FluidSystem::H2OIdx] = satVapPressure ;
479 beta[FluidSystem::N2Idx] = Henry ;
481 for (
int i = 0; i < numComponents; ++i)
482 priVars[conti00EqIdx + i] = xl[i]*beta[i];
490 bool onLeftBoundary_(
const GlobalPosition & globalPos)
const
491 {
return globalPos[0] < this->
gridGeometry().bBoxMin()[0] + eps_; }
496 bool onRightBoundary_(
const GlobalPosition & globalPos)
const
497 {
return globalPos[0] > this->
gridGeometry().bBoxMax()[0] - eps_; }
502 bool onLowerBoundary_(
const GlobalPosition & globalPos)
const
503 {
return globalPos[dimWorld-1] < this->
gridGeometry().bBoxMin()[dimWorld-1] + eps_; }
506 static constexpr Scalar eps_ = 1e-6;
507 Scalar percentOfEquil_ ;
510 std::string outputName_;
511 Scalar heatIntoSolid_;
513 Scalar SwPMInitial_ ;
514 Scalar SwFFInitial_ ;
518 Dune::ParameterTree inputParameters_;
519 Scalar x_[numPhases][numComponents] ;
523 Scalar massFluxInjectedPhase_ ;
525 std::shared_ptr<GridVariables> gridVariables_;
530 {
return inputParameters_; }
Defines a type tag and some properties for models using the box scheme.
Setting constant fluid properties via the input file.
Specification of a function relating volume specific interfacial area to capillary pressure and satur...
Specification of a function relating volume specific interfacial area to capillary pressure and satur...
Specification of a function relating volume specific interfacial area to capillary pressure and satur...
Specification of a function relating volume specific interfacial area to capillary pressure and satur...
Specification of a function relating volume specific interfacial area to capillary pressure and satur...
This material law takes a material law (interfacial area surface) defined for effective saturations a...
Enumeration of the formulations accepted by the MpNc model.
MpNcPressureFormulation
Enumerates the formulations which the MpNc model accepts.
Definition: pressureformulation.hh:36
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Property tag Indices
Definition: porousmediumflow/sequential/properties.hh:59
std::string molarDensity(int phaseIdx) noexcept
I/O name of molar density for multiphase systems.
Definition: name.hh:83
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Base class for all finite-volume problems.
Definition: common/fvproblem.hh:50
const GridGeometry & gridGeometry() const
The finite volume grid geometry.
Definition: common/fvproblem.hh:588
Property to specify the type of scalar values.
Definition: common/properties.hh:53
The DUNE grid type.
Definition: common/properties.hh:57
Property to specify the type of a problem which has to be solved.
Definition: common/properties.hh:69
The type of the spatial parameters object.
Definition: common/properties.hh:221
The type of the fluid system to use.
Definition: common/properties.hh:223
The type of the solid system to use.
Definition: common/properties.hh:227
Definition: common/properties.hh:312
static Scalar henry(Scalar temperature)
Henry coefficient for molecular nitrogen in liquid water.
Definition: h2o_n2.hh:48
A component which returns run time specified values for all fluid properties.
Definition: constant.hh:58
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: misciblemultiphasecomposition.hh:60
static void solve(FluidState &fluidState, ParameterCache ¶mCache, int knownPhaseIdx=0)
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: misciblemultiphasecomposition.hh:82
Policy for the H2O-N2 fluid system.
Definition: h2on2.hh:52
A two-phase fluid system with two components water Nitrogen for non-equilibrium models....
Definition: h2on2kinetic.hh:49
A solid phase consisting of a single inert solid component.
Definition: inertsolidphase.hh:41
Base class for all fully implicit porous media problems.
Definition: dumux/porousmediumflow/problem.hh:39
SpatialParams & spatialParams()
Returns the spatial parameters object.
Definition: dumux/porousmediumflow/problem.hh:146
Problem showcasing the capabilities of the kinetic model.
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:131
Dune::ParameterTree getInputParameters() const
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:529
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
Specifies which kind of boundary condition should be used for which equation on a given boundary segm...
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:230
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
Evaluates the boundary conditions for a Dirichlet boundary segment.
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:250
const GridVariables & gridVariables() const
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:201
EvaporationAtmosphereProblem(std::shared_ptr< const GridGeometry > gridGeometry)
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:178
NumEqVector neumann(const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const ElementFluxVariablesCache &elemFluxVarsCache, const SubControlVolumeFace &scvf) const
Evaluates the boundary conditions for a Neumann boundary segment.
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:258
NumEqVector source(const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolume &scv) const
Evaluates the source term for all balance equations within a given sub-controlvolume.
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:347
void setGridVariables(std::shared_ptr< GridVariables > gridVariables)
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:198
const std::string name() const
Returns the problem name.
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:214
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
Evaluates the initial value for a control volume.
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:331
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:71
std::tuple< MPNCNonequil > InheritsFrom
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:71
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:72
std::tuple< EvaporationAtmosphere, BoxModel > InheritsFrom
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:72
Dune::YaspGrid< 2, Dune::TensorProductCoordinates< GetPropType< TypeTag, Properties::Scalar >, 2 > > type
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:77
double type
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:101
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:107
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:117
GetPropType< TypeTag, GridGeometry > GridGeometry
Definition: test/porousmediumflow/mpnc/implicit/kinetic/problem.hh:116
Definition of the spatial parameters for the evaporation atmosphere Problem (using a "poor man's coup...
Definition: porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh:53
Declares all properties used in Dumux.
A fully implicit model for MpNc flow using vertex centered finite volumes.
Base class for all porous media problems.
Definition of the spatial parameters for the MaxwellStefan problem.