28 #ifndef EWOMS_WATER_AIR_PROBLEM_HH
29 #define EWOMS_WATER_AIR_PROBLEM_HH
34 #include <opm/material/fluidsystems/H2OAirFluidSystem.hpp>
35 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
36 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
37 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
38 #include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
39 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
40 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
41 #include <opm/material/heatconduction/Somerton.hpp>
42 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
43 #include <opm/common/Unused.hpp>
45 #include <dune/grid/yaspgrid.hh>
46 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
48 #include <dune/common/fvector.hh>
49 #include <dune/common/fmatrix.hh>
50 #include <dune/common/version.hh>
56 template <
class TypeTag>
61 namespace Properties {
71 SET_PROP(WaterAirBaseProblem, MaterialLaw)
75 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
76 typedef Opm::TwoPhaseMaterialTraits<Scalar,
77 FluidSystem::liquidPhaseIdx,
78 FluidSystem::gasPhaseIdx> Traits;
82 typedef Opm::RegularizedBrooksCorey<Traits> EffMaterialLaw;
87 typedef Opm::EffToAbsLaw<EffMaterialLaw> type;
91 SET_PROP(WaterAirBaseProblem, HeatConductionLaw)
95 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
99 typedef Opm::Somerton<FluidSystem, Scalar> type;
105 Opm::FluidSystems::H2OAir<typename
GET_PROP_TYPE(TypeTag, Scalar)>);
111 SET_INT_PROP(WaterAirBaseProblem, NumericDifferenceMethod, +1);
114 SET_BOOL_PROP(WaterAirBaseProblem, NewtonWriteConvergence, false);
117 SET_SCALAR_PROP(WaterAirBaseProblem, EndTime, 1.0 * 365 * 24 * 60 * 60);
126 SET_TAG_PROP(WaterAirBaseProblem, LinearSolverSplice, ParallelIstlLinearSolver);
128 Ewoms::Linear::SolverWrapperRestartedGMRes<TypeTag>);
130 Ewoms::Linear::PreconditionerWrapperILUn<TypeTag>);
131 SET_INT_PROP(WaterAirBaseProblem, PreconditionerOrder, 2);
164 template <
class TypeTag >
167 typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
173 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
176 numPhases = FluidSystem::numPhases,
179 temperatureIdx = Indices::temperatureIdx,
180 energyEqIdx = Indices::energyEqIdx,
183 H2OIdx = FluidSystem::H2OIdx,
184 AirIdx = FluidSystem::AirIdx,
187 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
188 gasPhaseIdx = FluidSystem::gasPhaseIdx,
191 conti0EqIdx = Indices::conti0EqIdx,
194 dim = GridView::dimension,
195 dimWorld = GridView::dimensionworld
198 static const bool enableEnergy =
GET_PROP_VALUE(TypeTag, EnableEnergy);
201 typedef typename
GET_PROP_TYPE(TypeTag, RateVector) RateVector;
202 typedef typename
GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
203 typedef typename
GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
204 typedef typename
GET_PROP_TYPE(TypeTag, Constraints) Constraints;
205 typedef typename
GET_PROP_TYPE(TypeTag, Simulator) Simulator;
207 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
208 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
209 typedef typename
GET_PROP_TYPE(TypeTag, HeatConductionLaw) HeatConductionLaw;
210 typedef typename
GET_PROP_TYPE(TypeTag, HeatConductionLawParams) HeatConductionLawParams;
212 typedef typename GridView::ctype CoordScalar;
213 typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
215 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
222 : ParentType(simulator)
230 ParentType::finishInit();
235 FluidSystem::init(275, 600, 100,
241 fineK_ = this->toDimMatrix_(1e-13);
242 coarseK_ = this->toDimMatrix_(1e-12);
246 coarsePorosity_ = 0.3;
249 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
250 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
251 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
252 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
255 fineMaterialParams_.setEntryPressure(1e4);
256 coarseMaterialParams_.setEntryPressure(1e4);
257 fineMaterialParams_.setLambda(2.0);
258 coarseMaterialParams_.setLambda(2.0);
260 fineMaterialParams_.finalize();
261 coarseMaterialParams_.finalize();
264 computeHeatCondParams_(fineHeatCondParams_, finePorosity_);
265 computeHeatCondParams_(coarseHeatCondParams_, coarsePorosity_);
278 std::ostringstream oss;
279 oss <<
"waterair_" << Model::name();
298 this->model().globalStorage(storage);
301 if (this->gridView().comm().rank() == 0) {
302 std::cout <<
"Storage: " << storage << std::endl << std::flush;
313 template <
class Context>
316 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
317 if (isFineMaterial_(pos))
325 template <
class Context>
326 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
328 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
329 if (isFineMaterial_(pos))
330 return finePorosity_;
332 return coarsePorosity_;
338 template <
class Context>
341 unsigned timeIdx)
const
343 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
344 if (isFineMaterial_(pos))
345 return fineMaterialParams_;
347 return coarseMaterialParams_;
355 template <
class Context>
357 unsigned spaceIdx OPM_UNUSED,
358 unsigned timeIdx OPM_UNUSED)
const
368 template <
class Context>
369 const HeatConductionLawParams&
372 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
373 if (isFineMaterial_(pos))
374 return fineHeatCondParams_;
375 return coarseHeatCondParams_;
393 template <
class Context>
395 const Context& context,
396 unsigned spaceIdx,
unsigned timeIdx)
const
398 const auto& pos = context.cvCenter(spaceIdx, timeIdx);
399 assert(onLeftBoundary_(pos) ||
400 onLowerBoundary_(pos) ||
401 onRightBoundary_(pos) ||
402 onUpperBoundary_(pos));
405 RateVector massRate(0.0);
406 massRate[conti0EqIdx + AirIdx] = -1e-3;
409 values.setMassRate(massRate);
412 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
413 initialFluidState_(fs, context, spaceIdx, timeIdx);
415 Scalar hl = fs.enthalpy(liquidPhaseIdx);
416 Scalar hg = fs.enthalpy(gasPhaseIdx);
417 values.setEnthalpyRate(values[conti0EqIdx + AirIdx] * hg +
418 values[conti0EqIdx + H2OIdx] * hl);
421 else if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
422 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
423 initialFluidState_(fs, context, spaceIdx, timeIdx);
426 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
446 template <
class Context>
448 const Context& context,
450 unsigned timeIdx)
const
452 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
453 initialFluidState_(fs, context, spaceIdx, timeIdx);
455 const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
456 values.assignMassConservative(fs, matParams,
true);
465 template <
class Context>
467 const Context& context OPM_UNUSED,
468 unsigned spaceIdx OPM_UNUSED,
469 unsigned timeIdx OPM_UNUSED)
const
475 bool onLeftBoundary_(
const GlobalPosition& pos)
const
476 {
return pos[0] < eps_; }
478 bool onRightBoundary_(
const GlobalPosition& pos)
const
479 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
481 bool onLowerBoundary_(
const GlobalPosition& pos)
const
482 {
return pos[1] < eps_; }
484 bool onUpperBoundary_(
const GlobalPosition& pos)
const
485 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
487 bool onInlet_(
const GlobalPosition& pos)
const
488 {
return onLowerBoundary_(pos) && (15.0 < pos[0]) && (pos[0] < 25.0); }
490 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const
491 {
return (20 < pos[0]) && (pos[0] < 30) && (pos[1] < 30); }
493 template <
class Context,
class Flu
idState>
494 void initialFluidState_(FluidState& fs,
495 const Context& context,
497 unsigned timeIdx)
const
499 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
501 Scalar densityW = 1000.0;
502 fs.setPressure(liquidPhaseIdx, 1e5 + (maxDepth_ - pos[1])*densityW*9.81);
503 fs.setSaturation(liquidPhaseIdx, 1.0);
504 fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
505 fs.setMoleFraction(liquidPhaseIdx, AirIdx, 0.0);
507 if (inHighTemperatureRegion_(pos))
508 fs.setTemperature(380);
510 fs.setTemperature(283.0 + (maxDepth_ - pos[1])*0.03);
513 fs.setSaturation(gasPhaseIdx, 0);
514 Scalar pc[numPhases];
515 const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
516 MaterialLaw::capillaryPressures(pc, matParams, fs);
517 fs.setPressure(gasPhaseIdx, fs.pressure(liquidPhaseIdx) + (pc[gasPhaseIdx] - pc[liquidPhaseIdx]));
519 typename FluidSystem::template ParameterCache<Scalar> paramCache;
520 typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> CFRP;
521 CFRP::solve(fs, paramCache, liquidPhaseIdx,
false,
true);
524 void computeHeatCondParams_(HeatConductionLawParams& params, Scalar poro)
526 Scalar lambdaGranite = 2.8;
529 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
530 fs.setTemperature(293.15);
531 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
532 fs.setPressure(phaseIdx, 1.0135e5);
535 typename FluidSystem::template ParameterCache<Scalar> paramCache;
536 paramCache.updateAll(fs);
537 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
538 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
539 fs.setDensity(phaseIdx, rho);
542 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
543 Scalar lambdaSaturated;
544 if (FluidSystem::isLiquid(phaseIdx)) {
546 FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
547 lambdaSaturated = std::pow(lambdaGranite, (1-poro)) + std::pow(lambdaFluid, poro);
550 lambdaSaturated = std::pow(lambdaGranite, (1-poro));
552 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
553 if (!FluidSystem::isLiquid(phaseIdx))
554 params.setVacuumLambda(lambdaSaturated);
558 bool isFineMaterial_(
const GlobalPosition& pos)
const
559 {
return pos[dim-1] > layerBottom_; }
565 Scalar finePorosity_;
566 Scalar coarsePorosity_;
568 MaterialLawParams fineMaterialParams_;
569 MaterialLawParams coarseMaterialParams_;
571 HeatConductionLawParams fineHeatCondParams_;
572 HeatConductionLawParams coarseHeatCondParams_;
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: waterairproblem.hh:447
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: waterairproblem.hh:228
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
Scalar heatCapacitySolid(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: waterairproblem.hh:356
std::string name() const
The problem name.
Definition: waterairproblem.hh:276
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:339
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
#define SET_TAG_PROP(EffTypeTagName, PropTagName, ValueTypeTagName)
Define a property containing a type tag.
Definition: propertysystem.hh:436
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: waterairproblem.hh:394
const HeatConductionLawParams & heatConductionParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:370
Provides all unmodified linear solvers from dune-istl.
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:326
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:314
Declares the properties required for the compositional multi-phase primary variable switching model...
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
#define SET_STRING_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant string value.
Definition: propertysystem.hh:416
Non-isothermal gas injection problem where a air is injected into a fully water saturated medium...
Definition: waterairproblem.hh:57
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: waterairproblem.hh:466
void endTimeStep()
Called by the simulator after each time integration.
Definition: waterairproblem.hh:289