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>
55template <
class TypeTag>
59namespace Opm::Properties {
66template<
class TypeTag>
67struct Grid<TypeTag, TTag::WaterAirBaseProblem> {
using type = Dune::YaspGrid<2>; };
70template<
class TypeTag>
74template<
class TypeTag>
75struct MaterialLaw<TypeTag, TTag::WaterAirBaseProblem>
78 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
79 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
80 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
81 FluidSystem::liquidPhaseIdx,
82 FluidSystem::gasPhaseIdx>;
86 using EffMaterialLaw = Opm::RegularizedBrooksCorey<Traits>;
91 using type = Opm::EffToAbsLaw<EffMaterialLaw>;
95template<
class TypeTag>
96struct ThermalConductionLaw<TypeTag, TTag::WaterAirBaseProblem>
99 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
100 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
104 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
108template<
class TypeTag>
109struct SolidEnergyLaw<TypeTag, TTag::WaterAirBaseProblem>
110{
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
114template<
class TypeTag>
115struct FluidSystem<TypeTag, TTag::WaterAirBaseProblem>
116{
using type = Opm::H2OAirFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
119template<
class TypeTag>
120struct EnableGravity<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr bool value =
true; };
123template<
class TypeTag>
124struct NumericDifferenceMethod<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr int value = +1; };
127template<
class TypeTag>
128struct NewtonWriteConvergence<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr bool value =
false; };
131template<
class TypeTag>
132struct EndTime<TypeTag, TTag::WaterAirBaseProblem>
134 using type = GetPropType<TypeTag, Scalar>;
135 static constexpr type value = 1.0 * 365 * 24 * 60 * 60;
139template<
class TypeTag>
140struct InitialTimeStepSize<TypeTag, TTag::WaterAirBaseProblem>
142 using type = GetPropType<TypeTag, Scalar>;
143 static constexpr type value = 250;
147template<
class TypeTag>
148struct GridFile<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr auto value =
"./data/waterair.dgf"; };
151template<
class TypeTag>
152struct LinearSolverSplice<TypeTag, TTag::WaterAirBaseProblem>
153{
using type = TTag::ParallelIstlLinearSolver; };
155template<
class TypeTag>
156struct LinearSolverWrapper<TypeTag, TTag::WaterAirBaseProblem>
157{
using type = Opm::Linear::SolverWrapperRestartedGMRes<TypeTag>; };
159template<
class TypeTag>
160struct PreconditionerWrapper<TypeTag, TTag::WaterAirBaseProblem>
161{
using type = Opm::Linear::PreconditionerWrapperILU<TypeTag>; };
162template<
class TypeTag>
163struct PreconditionerOrder<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr int value = 2; };
196template <
class TypeTag >
199 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
201 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
202 using GridView = GetPropType<TypeTag, Properties::GridView>;
205 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
206 using Indices = GetPropType<TypeTag, Properties::Indices>;
208 numPhases = FluidSystem::numPhases,
211 temperatureIdx = Indices::temperatureIdx,
212 energyEqIdx = Indices::energyEqIdx,
215 H2OIdx = FluidSystem::H2OIdx,
216 AirIdx = FluidSystem::AirIdx,
219 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
220 gasPhaseIdx = FluidSystem::gasPhaseIdx,
223 conti0EqIdx = Indices::conti0EqIdx,
226 dim = GridView::dimension,
227 dimWorld = GridView::dimensionworld
230 static const bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
232 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
233 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
234 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
235 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
236 using Constraints = GetPropType<TypeTag, Properties::Constraints>;
237 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
238 using Model = GetPropType<TypeTag, Properties::Model>;
239 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
240 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
241 using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
242 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
244 using CoordScalar =
typename GridView::ctype;
245 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
247 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
254 : ParentType(simulator)
262 ParentType::finishInit();
267 FluidSystem::init(275, 600, 100,
273 fineK_ = this->toDimMatrix_(1e-13);
274 coarseK_ = this->toDimMatrix_(1e-12);
278 coarsePorosity_ = 0.3;
281 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
282 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
283 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
284 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
287 fineMaterialParams_.setEntryPressure(1e4);
288 coarseMaterialParams_.setEntryPressure(1e4);
289 fineMaterialParams_.setLambda(2.0);
290 coarseMaterialParams_.setLambda(2.0);
292 fineMaterialParams_.finalize();
293 coarseMaterialParams_.finalize();
296 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
297 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
300 solidEnergyLawParams_.setSolidHeatCapacity(790.0
302 solidEnergyLawParams_.finalize();
315 std::ostringstream oss;
316 oss <<
"waterair_" << Model::name();
317 if (getPropValue<TypeTag, Properties::EnableEnergy>())
335 this->model().globalStorage(storage);
338 if (this->gridView().comm().rank() == 0) {
339 std::cout <<
"Storage: " << storage << std::endl << std::flush;
350 template <
class Context>
353 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
354 if (isFineMaterial_(pos))
362 template <
class Context>
363 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
365 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
366 if (isFineMaterial_(pos))
367 return finePorosity_;
369 return coarsePorosity_;
375 template <
class Context>
378 unsigned timeIdx)
const
380 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
381 if (isFineMaterial_(pos))
382 return fineMaterialParams_;
384 return coarseMaterialParams_;
392 template <
class Context>
393 const SolidEnergyLawParams&
397 {
return solidEnergyLawParams_; }
402 template <
class Context>
403 const ThermalConductionLawParams&
406 unsigned timeIdx)
const
408 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
409 if (isFineMaterial_(pos))
410 return fineThermalCondParams_;
411 return coarseThermalCondParams_;
429 template <
class Context>
431 const Context& context,
432 unsigned spaceIdx,
unsigned timeIdx)
const
434 const auto& pos = context.cvCenter(spaceIdx, timeIdx);
435 assert(onLeftBoundary_(pos) ||
436 onLowerBoundary_(pos) ||
437 onRightBoundary_(pos) ||
438 onUpperBoundary_(pos));
441 RateVector massRate(0.0);
442 massRate[conti0EqIdx + AirIdx] = -1e-3;
445 values.setMassRate(massRate);
448 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
449 initialFluidState_(fs, context, spaceIdx, timeIdx);
451 Scalar hl = fs.enthalpy(liquidPhaseIdx);
452 Scalar hg = fs.enthalpy(gasPhaseIdx);
453 values.setEnthalpyRate(values[conti0EqIdx + AirIdx] * hg +
454 values[conti0EqIdx + H2OIdx] * hl);
457 else if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
458 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
459 initialFluidState_(fs, context, spaceIdx, timeIdx);
462 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
482 template <
class Context>
484 const Context& context,
486 unsigned timeIdx)
const
488 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
489 initialFluidState_(fs, context, spaceIdx, timeIdx);
492 values.assignMassConservative(fs, matParams,
true);
501 template <
class Context>
511 bool onLeftBoundary_(
const GlobalPosition& pos)
const
512 {
return pos[0] < eps_; }
514 bool onRightBoundary_(
const GlobalPosition& pos)
const
515 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
517 bool onLowerBoundary_(
const GlobalPosition& pos)
const
518 {
return pos[1] < eps_; }
520 bool onUpperBoundary_(
const GlobalPosition& pos)
const
521 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
523 bool onInlet_(
const GlobalPosition& pos)
const
524 {
return onLowerBoundary_(pos) && (15.0 < pos[0]) && (pos[0] < 25.0); }
526 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const
527 {
return (20 < pos[0]) && (pos[0] < 30) && (pos[1] < 30); }
529 template <
class Context,
class Flu
idState>
530 void initialFluidState_(FluidState& fs,
531 const Context& context,
533 unsigned timeIdx)
const
535 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
537 Scalar densityW = 1000.0;
538 fs.setPressure(liquidPhaseIdx, 1e5 + (maxDepth_ - pos[1])*densityW*9.81);
539 fs.setSaturation(liquidPhaseIdx, 1.0);
540 fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
541 fs.setMoleFraction(liquidPhaseIdx, AirIdx, 0.0);
543 if (inHighTemperatureRegion_(pos))
544 fs.setTemperature(380);
546 fs.setTemperature(283.0 + (maxDepth_ - pos[1])*0.03);
549 fs.setSaturation(gasPhaseIdx, 0);
550 Scalar pc[numPhases];
552 MaterialLaw::capillaryPressures(pc, matParams, fs);
553 fs.setPressure(gasPhaseIdx, fs.pressure(liquidPhaseIdx) + (pc[gasPhaseIdx] - pc[liquidPhaseIdx]));
555 typename FluidSystem::template ParameterCache<Scalar> paramCache;
556 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
557 CFRP::solve(fs, paramCache, liquidPhaseIdx,
true,
true);
560 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
562 Scalar lambdaGranite = 2.8;
565 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
566 fs.setTemperature(293.15);
567 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
568 fs.setPressure(phaseIdx, 1.0135e5);
571 typename FluidSystem::template ParameterCache<Scalar> paramCache;
572 paramCache.updateAll(fs);
573 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
574 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
575 fs.setDensity(phaseIdx, rho);
578 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
579 Scalar lambdaSaturated;
580 if (FluidSystem::isLiquid(phaseIdx)) {
582 FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
583 lambdaSaturated = std::pow(lambdaGranite, (1-poro)) + std::pow(lambdaFluid, poro);
586 lambdaSaturated = std::pow(lambdaGranite, (1-poro));
588 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
589 if (!FluidSystem::isLiquid(phaseIdx))
590 params.setVacuumLambda(lambdaSaturated);
594 bool isFineMaterial_(
const GlobalPosition& pos)
const
595 {
return pos[dim-1] > layerBottom_; }
601 Scalar finePorosity_;
602 Scalar coarsePorosity_;
604 MaterialLawParams fineMaterialParams_;
605 MaterialLawParams coarseMaterialParams_;
607 ThermalConductionLawParams fineThermalCondParams_;
608 ThermalConductionLawParams coarseThermalCondParams_;
609 SolidEnergyLawParams solidEnergyLawParams_;
Non-isothermal gas injection problem where a air is injected into a fully water saturated medium.
Definition: waterairproblem.hh:198
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:351
void finishInit()
Definition: waterairproblem.hh:260
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:483
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:430
std::string name() const
Definition: waterairproblem.hh:313
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:363
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: waterairproblem.hh:502
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:404
void endTimeStep()
Definition: waterairproblem.hh:326
WaterAirProblem(Simulator &simulator)
Definition: waterairproblem.hh:253
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:376
const SolidEnergyLawParams & solidEnergyLawParams(const Context &, unsigned, unsigned) const
Return the parameters for the energy storage law of the rock.
Definition: waterairproblem.hh:394
Definition: waterairproblem.hh:62