28#ifndef EWOMS_WATER_AIR_PROBLEM_HH
29#define EWOMS_WATER_AIR_PROBLEM_HH
31#include <dune/common/fmatrix.hh>
32#include <dune/common/fvector.hh>
34#include <dune/grid/yaspgrid.hh>
35#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
37#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
38#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
39#include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
40#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
41#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
42#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
43#include <opm/material/fluidstates/CompositionalFluidState.hpp>
44#include <opm/material/fluidsystems/H2OAirFluidSystem.hpp>
45#include <opm/material/thermal/ConstantSolidHeatCapLaw.hpp>
46#include <opm/material/thermal/SomertonThermalConductionLaw.hpp>
60template <
class TypeTag>
71template<
class TypeTag>
72struct Grid<TypeTag, TTag::WaterAirBaseProblem> {
using type = Dune::YaspGrid<2>; };
75template<
class TypeTag>
79template<
class TypeTag>
85 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
86 FluidSystem::liquidPhaseIdx,
87 FluidSystem::gasPhaseIdx>;
91 using EffMaterialLaw = Opm::RegularizedBrooksCorey<Traits>;
96 using type = Opm::EffToAbsLaw<EffMaterialLaw>;
100template<
class TypeTag>
109 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
113template<
class TypeTag>
115{
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
119template<
class TypeTag>
121{
using type = Opm::H2OAirFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
124template<
class TypeTag>
128template<
class TypeTag>
132template<
class TypeTag>
167template <
class TypeTag >
179 numPhases = FluidSystem::numPhases,
182 temperatureIdx = Indices::temperatureIdx,
183 energyEqIdx = Indices::energyEqIdx,
186 H2OIdx = FluidSystem::H2OIdx,
187 AirIdx = FluidSystem::AirIdx,
190 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
191 gasPhaseIdx = FluidSystem::gasPhaseIdx,
194 conti0EqIdx = Indices::conti0EqIdx,
197 dim = GridView::dimension,
198 dimWorld = GridView::dimensionworld
201 static const bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
215 using CoordScalar =
typename GridView::ctype;
216 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
218 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
225 : ParentType(simulator)
233 ParentType::finishInit();
238 FluidSystem::init(275, 600, 100,
244 fineK_ = this->toDimMatrix_(1e-13);
245 coarseK_ = this->toDimMatrix_(1e-12);
249 coarsePorosity_ = 0.3;
252 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
253 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
254 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
255 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
258 fineMaterialParams_.setEntryPressure(1e4);
259 coarseMaterialParams_.setEntryPressure(1e4);
260 fineMaterialParams_.setLambda(2.0);
261 coarseMaterialParams_.setLambda(2.0);
263 fineMaterialParams_.finalize();
264 coarseMaterialParams_.finalize();
267 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
268 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
271 solidEnergyLawParams_.setSolidHeatCapacity(790.0
273 solidEnergyLawParams_.finalize();
285 ParentType::registerParameters();
287 Parameters::SetDefault<Parameters::GridFile>(
"./data/waterair.dgf");
290 Parameters::SetDefault<Parameters::NumericDifferenceMethod>(+1);
292 Parameters::SetDefault<Parameters::EndTime<Scalar>>(1.0 * 365 * 24 * 60 * 60);
293 Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(250.0);
294 Parameters::SetDefault<Parameters::EnableGravity>(
true);
295 Parameters::SetDefault<Parameters::PreconditionerOrder>(2);
303 std::ostringstream oss;
304 oss <<
"waterair_" << Model::name();
305 if (getPropValue<TypeTag, Properties::EnableEnergy>())
323 this->model().globalStorage(storage);
326 if (this->gridView().comm().rank() == 0) {
327 std::cout <<
"Storage: " << storage << std::endl << std::flush;
338 template <
class Context>
341 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
342 if (isFineMaterial_(pos))
350 template <
class Context>
351 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
353 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
354 if (isFineMaterial_(pos))
355 return finePorosity_;
357 return coarsePorosity_;
363 template <
class Context>
366 unsigned timeIdx)
const
368 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
369 if (isFineMaterial_(pos))
370 return fineMaterialParams_;
372 return coarseMaterialParams_;
380 template <
class Context>
381 const SolidEnergyLawParams&
385 {
return solidEnergyLawParams_; }
390 template <
class Context>
391 const ThermalConductionLawParams&
394 unsigned timeIdx)
const
396 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
397 if (isFineMaterial_(pos))
398 return fineThermalCondParams_;
399 return coarseThermalCondParams_;
417 template <
class Context>
419 const Context& context,
420 unsigned spaceIdx,
unsigned timeIdx)
const
422 const auto& pos = context.cvCenter(spaceIdx, timeIdx);
423 assert(onLeftBoundary_(pos) ||
424 onLowerBoundary_(pos) ||
425 onRightBoundary_(pos) ||
426 onUpperBoundary_(pos));
429 RateVector massRate(0.0);
430 massRate[conti0EqIdx + AirIdx] = -1e-3;
433 values.setMassRate(massRate);
436 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
437 initialFluidState_(fs, context, spaceIdx, timeIdx);
439 Scalar hl = fs.enthalpy(liquidPhaseIdx);
440 Scalar hg = fs.enthalpy(gasPhaseIdx);
441 values.setEnthalpyRate(values[conti0EqIdx + AirIdx] * hg +
442 values[conti0EqIdx + H2OIdx] * hl);
445 else if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
446 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
447 initialFluidState_(fs, context, spaceIdx, timeIdx);
450 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
470 template <
class Context>
472 const Context& context,
474 unsigned timeIdx)
const
476 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
477 initialFluidState_(fs, context, spaceIdx, timeIdx);
480 values.assignMassConservative(fs, matParams,
true);
489 template <
class Context>
499 bool onLeftBoundary_(
const GlobalPosition& pos)
const
500 {
return pos[0] < eps_; }
502 bool onRightBoundary_(
const GlobalPosition& pos)
const
503 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
505 bool onLowerBoundary_(
const GlobalPosition& pos)
const
506 {
return pos[1] < eps_; }
508 bool onUpperBoundary_(
const GlobalPosition& pos)
const
509 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
511 bool onInlet_(
const GlobalPosition& pos)
const
512 {
return onLowerBoundary_(pos) && (15.0 < pos[0]) && (pos[0] < 25.0); }
514 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const
515 {
return (20 < pos[0]) && (pos[0] < 30) && (pos[1] < 30); }
517 template <
class Context,
class Flu
idState>
518 void initialFluidState_(FluidState& fs,
519 const Context& context,
521 unsigned timeIdx)
const
523 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
525 Scalar densityW = 1000.0;
526 fs.setPressure(liquidPhaseIdx, 1e5 + (maxDepth_ - pos[1])*densityW*9.81);
527 fs.setSaturation(liquidPhaseIdx, 1.0);
528 fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
529 fs.setMoleFraction(liquidPhaseIdx, AirIdx, 0.0);
531 if (inHighTemperatureRegion_(pos))
532 fs.setTemperature(380);
534 fs.setTemperature(283.0 + (maxDepth_ - pos[1])*0.03);
537 fs.setSaturation(gasPhaseIdx, 0);
538 Scalar pc[numPhases];
540 MaterialLaw::capillaryPressures(pc, matParams, fs);
541 fs.setPressure(gasPhaseIdx, fs.pressure(liquidPhaseIdx) + (pc[gasPhaseIdx] - pc[liquidPhaseIdx]));
543 typename FluidSystem::template ParameterCache<Scalar> paramCache;
544 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
545 CFRP::solve(fs, paramCache, liquidPhaseIdx,
true,
true);
548 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
550 Scalar lambdaGranite = 2.8;
553 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
554 fs.setTemperature(293.15);
555 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
556 fs.setPressure(phaseIdx, 1.0135e5);
559 typename FluidSystem::template ParameterCache<Scalar> paramCache;
560 paramCache.updateAll(fs);
561 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
562 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
563 fs.setDensity(phaseIdx, rho);
566 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
567 Scalar lambdaSaturated;
568 if (FluidSystem::isLiquid(phaseIdx)) {
570 FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
571 lambdaSaturated = std::pow(lambdaGranite, (1-poro)) + std::pow(lambdaFluid, poro);
574 lambdaSaturated = std::pow(lambdaGranite, (1-poro));
576 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
577 if (!FluidSystem::isLiquid(phaseIdx))
578 params.setVacuumLambda(lambdaSaturated);
582 bool isFineMaterial_(
const GlobalPosition& pos)
const
583 {
return pos[dim-1] > layerBottom_; }
589 Scalar finePorosity_;
590 Scalar coarsePorosity_;
592 MaterialLawParams fineMaterialParams_;
593 MaterialLawParams coarseMaterialParams_;
595 ThermalConductionLawParams fineThermalCondParams_;
596 ThermalConductionLawParams coarseThermalCondParams_;
597 SolidEnergyLawParams solidEnergyLawParams_;
Definition: istlpreconditionerwrappers.hh:153
Solver wrapper for the restarted GMRES solver of dune-istl.
Definition: istlsolverwrappers.hh:115
Non-isothermal gas injection problem where a air is injected into a fully water saturated medium.
Definition: waterairproblem.hh:169
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:339
static void registerParameters()
Definition: waterairproblem.hh:283
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: waterairproblem.hh:231
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: waterairproblem.hh:471
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: waterairproblem.hh:418
std::string name() const
The problem name.
Definition: waterairproblem.hh:301
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:351
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: waterairproblem.hh:490
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:392
void endTimeStep()
Called by the simulator after each time integration.
Definition: waterairproblem.hh:314
WaterAirProblem(Simulator &simulator)
Definition: waterairproblem.hh:224
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:364
const SolidEnergyLawParams & solidEnergyLawParams(const Context &, unsigned, unsigned) const
Return the parameters for the energy storage law of the rock.
Definition: waterairproblem.hh:382
Defines the common parameters for the porous medium multi-phase models.
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
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:235
Declares the properties required for the compositional multi-phase primary variable switching model.
Opm::H2OAirFluidSystem< GetPropType< TypeTag, Properties::Scalar > > type
Definition: waterairproblem.hh:121
The fluid systems including the information about the phases.
Definition: multiphasebaseproperties.hh:69
Dune::YaspGrid< 2 > type
Definition: waterairproblem.hh:72
The type of the DUNE grid.
Definition: basicproperties.hh:100
Definition: fvbaseproperties.hh:53
Definition: linalgproperties.hh:57
Opm::EffToAbsLaw< EffMaterialLaw > type
Definition: waterairproblem.hh:96
The material law which ought to be used (extracted from the spatial parameters)
Definition: multiphasebaseproperties.hh:51
the preconditioner used by the linear solver
Definition: linalgproperties.hh:42
The type of the problem.
Definition: fvbaseproperties.hh:81
Opm::ConstantSolidHeatCapLaw< GetPropType< TypeTag, Properties::Scalar > > type
Definition: waterairproblem.hh:115
The material law for the energy stored in the solid matrix.
Definition: multiphasebaseproperties.hh:57
Definition: parallelistlbackend.hh:40
Definition: waterairproblem.hh:67
Opm::SomertonThermalConductionLaw< FluidSystem, Scalar > type
Definition: waterairproblem.hh:109
The material law for thermal conduction.
Definition: multiphasebaseproperties.hh:63