28#ifndef EWOMS_WATER_AIR_PROBLEM_HH
29#define EWOMS_WATER_AIR_PROBLEM_HH
31#include <opm/models/pvs/pvsproperties.hh>
32#include <opm/simulators/linalg/parallelistlbackend.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/thermal/ConstantSolidHeatCapLaw.hpp>
42#include <opm/material/thermal/SomertonThermalConductionLaw.hpp>
43#include <opm/material/constraintsolvers/ComputeFromReferencePhase.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>
56template <
class TypeTag>
60namespace Opm::Properties {
67template<
class TypeTag>
68struct Grid<TypeTag, TTag::WaterAirBaseProblem> {
using type = Dune::YaspGrid<2>; };
71template<
class TypeTag>
75template<
class TypeTag>
76struct MaterialLaw<TypeTag, TTag::WaterAirBaseProblem>
79 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
80 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
81 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
82 FluidSystem::liquidPhaseIdx,
83 FluidSystem::gasPhaseIdx>;
87 using EffMaterialLaw = Opm::RegularizedBrooksCorey<Traits>;
92 using type = Opm::EffToAbsLaw<EffMaterialLaw>;
96template<
class TypeTag>
97struct ThermalConductionLaw<TypeTag, TTag::WaterAirBaseProblem>
100 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
101 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
105 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
109template<
class TypeTag>
110struct SolidEnergyLaw<TypeTag, TTag::WaterAirBaseProblem>
111{
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
115template<
class TypeTag>
116struct FluidSystem<TypeTag, TTag::WaterAirBaseProblem>
117{
using type = Opm::H2OAirFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
120template<
class TypeTag>
121struct EnableGravity<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr bool value =
true; };
124template<
class TypeTag>
125struct NumericDifferenceMethod<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr int value = +1; };
128template<
class TypeTag>
129struct NewtonWriteConvergence<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr bool value =
false; };
132template<
class TypeTag>
133struct EndTime<TypeTag, TTag::WaterAirBaseProblem>
135 using type = GetPropType<TypeTag, Scalar>;
136 static constexpr type value = 1.0 * 365 * 24 * 60 * 60;
140template<
class TypeTag>
141struct InitialTimeStepSize<TypeTag, TTag::WaterAirBaseProblem>
143 using type = GetPropType<TypeTag, Scalar>;
144 static constexpr type value = 250;
148template<
class TypeTag>
149struct GridFile<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr auto value =
"./data/waterair.dgf"; };
152template<
class TypeTag>
153struct LinearSolverSplice<TypeTag, TTag::WaterAirBaseProblem>
154{
using type = TTag::ParallelIstlLinearSolver; };
156template<
class TypeTag>
157struct LinearSolverWrapper<TypeTag, TTag::WaterAirBaseProblem>
158{
using type = Opm::Linear::SolverWrapperRestartedGMRes<TypeTag>; };
160#if DUNE_VERSION_NEWER(DUNE_ISTL, 2,7)
161template<
class TypeTag>
162struct PreconditionerWrapper<TypeTag, TTag::WaterAirBaseProblem>
163{
using type = Opm::Linear::PreconditionerWrapperILU<TypeTag>; };
165template<
class TypeTag>
166struct PreconditionerWrapper<TypeTag, TTag::WaterAirBaseProblem>
167{
using type = Opm::Linear::PreconditionerWrapperILUn<TypeTag>; };
169template<
class TypeTag>
170struct PreconditionerOrder<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr int value = 2; };
203template <
class TypeTag >
206 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
208 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
209 using GridView = GetPropType<TypeTag, Properties::GridView>;
212 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
213 using Indices = GetPropType<TypeTag, Properties::Indices>;
215 numPhases = FluidSystem::numPhases,
218 temperatureIdx = Indices::temperatureIdx,
219 energyEqIdx = Indices::energyEqIdx,
222 H2OIdx = FluidSystem::H2OIdx,
223 AirIdx = FluidSystem::AirIdx,
226 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
227 gasPhaseIdx = FluidSystem::gasPhaseIdx,
230 conti0EqIdx = Indices::conti0EqIdx,
233 dim = GridView::dimension,
234 dimWorld = GridView::dimensionworld
237 static const bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
239 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
240 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
241 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
242 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
243 using Constraints = GetPropType<TypeTag, Properties::Constraints>;
244 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
245 using Model = GetPropType<TypeTag, Properties::Model>;
246 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
247 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
248 using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
249 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
251 using CoordScalar =
typename GridView::ctype;
252 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
254 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
261 : ParentType(simulator)
269 ParentType::finishInit();
274 FluidSystem::init(275, 600, 100,
280 fineK_ = this->toDimMatrix_(1e-13);
281 coarseK_ = this->toDimMatrix_(1e-12);
285 coarsePorosity_ = 0.3;
288 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
289 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
290 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
291 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
294 fineMaterialParams_.setEntryPressure(1e4);
295 coarseMaterialParams_.setEntryPressure(1e4);
296 fineMaterialParams_.setLambda(2.0);
297 coarseMaterialParams_.setLambda(2.0);
299 fineMaterialParams_.finalize();
300 coarseMaterialParams_.finalize();
303 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
304 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
307 solidEnergyLawParams_.setSolidHeatCapacity(790.0
309 solidEnergyLawParams_.finalize();
322 std::ostringstream oss;
323 oss <<
"waterair_" << Model::name();
324 if (getPropValue<TypeTag, Properties::EnableEnergy>())
342 this->model().globalStorage(storage);
345 if (this->gridView().comm().rank() == 0) {
346 std::cout <<
"Storage: " << storage << std::endl << std::flush;
357 template <
class Context>
360 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
361 if (isFineMaterial_(pos))
369 template <
class Context>
370 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
372 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
373 if (isFineMaterial_(pos))
374 return finePorosity_;
376 return coarsePorosity_;
382 template <
class Context>
385 unsigned timeIdx)
const
387 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
388 if (isFineMaterial_(pos))
389 return fineMaterialParams_;
391 return coarseMaterialParams_;
399 template <
class Context>
400 const SolidEnergyLawParams&
404 {
return solidEnergyLawParams_; }
409 template <
class Context>
410 const ThermalConductionLawParams&
413 unsigned timeIdx)
const
415 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
416 if (isFineMaterial_(pos))
417 return fineThermalCondParams_;
418 return coarseThermalCondParams_;
436 template <
class Context>
438 const Context& context,
439 unsigned spaceIdx,
unsigned timeIdx)
const
441 const auto& pos = context.cvCenter(spaceIdx, timeIdx);
442 assert(onLeftBoundary_(pos) ||
443 onLowerBoundary_(pos) ||
444 onRightBoundary_(pos) ||
445 onUpperBoundary_(pos));
448 RateVector massRate(0.0);
449 massRate[conti0EqIdx + AirIdx] = -1e-3;
452 values.setMassRate(massRate);
455 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
456 initialFluidState_(fs, context, spaceIdx, timeIdx);
458 Scalar hl = fs.enthalpy(liquidPhaseIdx);
459 Scalar hg = fs.enthalpy(gasPhaseIdx);
460 values.setEnthalpyRate(values[conti0EqIdx + AirIdx] * hg +
461 values[conti0EqIdx + H2OIdx] * hl);
464 else if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
465 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
466 initialFluidState_(fs, context, spaceIdx, timeIdx);
469 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
489 template <
class Context>
491 const Context& context,
493 unsigned timeIdx)
const
495 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
496 initialFluidState_(fs, context, spaceIdx, timeIdx);
499 values.assignMassConservative(fs, matParams,
true);
508 template <
class Context>
518 bool onLeftBoundary_(
const GlobalPosition& pos)
const
519 {
return pos[0] < eps_; }
521 bool onRightBoundary_(
const GlobalPosition& pos)
const
522 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
524 bool onLowerBoundary_(
const GlobalPosition& pos)
const
525 {
return pos[1] < eps_; }
527 bool onUpperBoundary_(
const GlobalPosition& pos)
const
528 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
530 bool onInlet_(
const GlobalPosition& pos)
const
531 {
return onLowerBoundary_(pos) && (15.0 < pos[0]) && (pos[0] < 25.0); }
533 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const
534 {
return (20 < pos[0]) && (pos[0] < 30) && (pos[1] < 30); }
536 template <
class Context,
class Flu
idState>
537 void initialFluidState_(FluidState& fs,
538 const Context& context,
540 unsigned timeIdx)
const
542 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
544 Scalar densityW = 1000.0;
545 fs.setPressure(liquidPhaseIdx, 1e5 + (maxDepth_ - pos[1])*densityW*9.81);
546 fs.setSaturation(liquidPhaseIdx, 1.0);
547 fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
548 fs.setMoleFraction(liquidPhaseIdx, AirIdx, 0.0);
550 if (inHighTemperatureRegion_(pos))
551 fs.setTemperature(380);
553 fs.setTemperature(283.0 + (maxDepth_ - pos[1])*0.03);
556 fs.setSaturation(gasPhaseIdx, 0);
557 Scalar pc[numPhases];
559 MaterialLaw::capillaryPressures(pc, matParams, fs);
560 fs.setPressure(gasPhaseIdx, fs.pressure(liquidPhaseIdx) + (pc[gasPhaseIdx] - pc[liquidPhaseIdx]));
562 typename FluidSystem::template ParameterCache<Scalar> paramCache;
563 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
564 CFRP::solve(fs, paramCache, liquidPhaseIdx,
true,
true);
567 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
569 Scalar lambdaGranite = 2.8;
572 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
573 fs.setTemperature(293.15);
574 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
575 fs.setPressure(phaseIdx, 1.0135e5);
578 typename FluidSystem::template ParameterCache<Scalar> paramCache;
579 paramCache.updateAll(fs);
580 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
581 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
582 fs.setDensity(phaseIdx, rho);
585 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
586 Scalar lambdaSaturated;
587 if (FluidSystem::isLiquid(phaseIdx)) {
589 FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
590 lambdaSaturated = std::pow(lambdaGranite, (1-poro)) + std::pow(lambdaFluid, poro);
593 lambdaSaturated = std::pow(lambdaGranite, (1-poro));
595 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
596 if (!FluidSystem::isLiquid(phaseIdx))
597 params.setVacuumLambda(lambdaSaturated);
601 bool isFineMaterial_(
const GlobalPosition& pos)
const
602 {
return pos[dim-1] > layerBottom_; }
608 Scalar finePorosity_;
609 Scalar coarsePorosity_;
611 MaterialLawParams fineMaterialParams_;
612 MaterialLawParams coarseMaterialParams_;
614 ThermalConductionLawParams fineThermalCondParams_;
615 ThermalConductionLawParams coarseThermalCondParams_;
616 SolidEnergyLawParams solidEnergyLawParams_;
Non-isothermal gas injection problem where a air is injected into a fully water saturated medium.
Definition: waterairproblem.hh:205
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:358
void finishInit()
Definition: waterairproblem.hh:267
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:490
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:437
std::string name() const
Definition: waterairproblem.hh:320
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:370
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: waterairproblem.hh:509
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:411
void endTimeStep()
Definition: waterairproblem.hh:333
WaterAirProblem(Simulator &simulator)
Definition: waterairproblem.hh:260
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:383
const SolidEnergyLawParams & solidEnergyLawParams(const Context &, unsigned, unsigned) const
Return the parameters for the energy storage law of the rock.
Definition: waterairproblem.hh:401
Definition: waterairproblem.hh:63