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>
44 #include <opm/material/common/Unused.hpp>
46 #include <dune/grid/yaspgrid.hh>
47 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
49 #include <dune/common/fvector.hh>
50 #include <dune/common/fmatrix.hh>
51 #include <dune/common/version.hh>
57 template <
class TypeTag>
58 class WaterAirProblem;
61 namespace Opm::Properties {
68 template<
class TypeTag>
69 struct Grid<TypeTag, TTag::WaterAirBaseProblem> {
using type = Dune::YaspGrid<2>; };
72 template<
class TypeTag>
76 template<
class TypeTag>
77 struct MaterialLaw<TypeTag, TTag::WaterAirBaseProblem>
80 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
81 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
82 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
83 FluidSystem::liquidPhaseIdx,
84 FluidSystem::gasPhaseIdx>;
88 using EffMaterialLaw = Opm::RegularizedBrooksCorey<Traits>;
93 using type = Opm::EffToAbsLaw<EffMaterialLaw>;
97 template<
class TypeTag>
98 struct ThermalConductionLaw<TypeTag, TTag::WaterAirBaseProblem>
101 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
102 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
106 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
110 template<
class TypeTag>
111 struct SolidEnergyLaw<TypeTag, TTag::WaterAirBaseProblem>
112 {
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
116 template<
class TypeTag>
117 struct FluidSystem<TypeTag, TTag::WaterAirBaseProblem>
118 {
using type = Opm::H2OAirFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
121 template<
class TypeTag>
122 struct EnableGravity<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr
bool value =
true; };
125 template<
class TypeTag>
126 struct NumericDifferenceMethod<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr
int value = +1; };
129 template<
class TypeTag>
130 struct NewtonWriteConvergence<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr
bool value =
false; };
133 template<
class TypeTag>
134 struct EndTime<TypeTag, TTag::WaterAirBaseProblem>
136 using type = GetPropType<TypeTag, Scalar>;
137 static constexpr type value = 1.0 * 365 * 24 * 60 * 60;
141 template<
class TypeTag>
142 struct InitialTimeStepSize<TypeTag, TTag::WaterAirBaseProblem>
144 using type = GetPropType<TypeTag, Scalar>;
145 static constexpr type value = 250;
149 template<
class TypeTag>
150 struct GridFile<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr
auto value =
"./data/waterair.dgf"; };
153 template<
class TypeTag>
154 struct LinearSolverSplice<TypeTag, TTag::WaterAirBaseProblem>
155 {
using type = TTag::ParallelIstlLinearSolver; };
157 template<
class TypeTag>
158 struct LinearSolverWrapper<TypeTag, TTag::WaterAirBaseProblem>
159 {
using type = Opm::Linear::SolverWrapperRestartedGMRes<TypeTag>; };
161 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2,7)
162 template<
class TypeTag>
163 struct PreconditionerWrapper<TypeTag, TTag::WaterAirBaseProblem>
164 {
using type = Opm::Linear::PreconditionerWrapperILU<TypeTag>; };
166 template<
class TypeTag>
167 struct PreconditionerWrapper<TypeTag, TTag::WaterAirBaseProblem>
168 {
using type = Opm::Linear::PreconditionerWrapperILUn<TypeTag>; };
170 template<
class TypeTag>
171 struct PreconditionerOrder<TypeTag, TTag::WaterAirBaseProblem> {
static constexpr
int value = 2; };
204 template <
class TypeTag >
207 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
209 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
210 using GridView = GetPropType<TypeTag, Properties::GridView>;
213 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
214 using Indices = GetPropType<TypeTag, Properties::Indices>;
216 numPhases = FluidSystem::numPhases,
219 temperatureIdx = Indices::temperatureIdx,
220 energyEqIdx = Indices::energyEqIdx,
223 H2OIdx = FluidSystem::H2OIdx,
224 AirIdx = FluidSystem::AirIdx,
227 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
228 gasPhaseIdx = FluidSystem::gasPhaseIdx,
231 conti0EqIdx = Indices::conti0EqIdx,
234 dim = GridView::dimension,
235 dimWorld = GridView::dimensionworld
238 static const bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
240 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
241 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
242 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
243 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
244 using Constraints = GetPropType<TypeTag, Properties::Constraints>;
245 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
246 using Model = GetPropType<TypeTag, Properties::Model>;
247 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
248 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
249 using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
250 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
252 using CoordScalar =
typename GridView::ctype;
253 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
255 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
262 : ParentType(simulator)
270 ParentType::finishInit();
275 FluidSystem::init(275, 600, 100,
281 fineK_ = this->toDimMatrix_(1e-13);
282 coarseK_ = this->toDimMatrix_(1e-12);
286 coarsePorosity_ = 0.3;
289 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
290 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
291 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
292 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
295 fineMaterialParams_.setEntryPressure(1e4);
296 coarseMaterialParams_.setEntryPressure(1e4);
297 fineMaterialParams_.setLambda(2.0);
298 coarseMaterialParams_.setLambda(2.0);
300 fineMaterialParams_.finalize();
301 coarseMaterialParams_.finalize();
304 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
305 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
308 solidEnergyLawParams_.setSolidHeatCapacity(790.0
310 solidEnergyLawParams_.finalize();
323 std::ostringstream oss;
324 oss <<
"waterair_" << Model::name();
325 if (getPropValue<TypeTag, Properties::EnableEnergy>())
343 this->model().globalStorage(storage);
346 if (this->gridView().comm().rank() == 0) {
347 std::cout <<
"Storage: " << storage << std::endl << std::flush;
358 template <
class Context>
361 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
362 if (isFineMaterial_(pos))
370 template <
class Context>
371 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
373 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
374 if (isFineMaterial_(pos))
375 return finePorosity_;
377 return coarsePorosity_;
383 template <
class Context>
386 unsigned timeIdx)
const
388 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
389 if (isFineMaterial_(pos))
390 return fineMaterialParams_;
392 return coarseMaterialParams_;
400 template <
class Context>
401 const SolidEnergyLawParams&
403 unsigned spaceIdx OPM_UNUSED,
404 unsigned timeIdx OPM_UNUSED)
const
405 {
return solidEnergyLawParams_; }
410 template <
class Context>
411 const ThermalConductionLawParams&
414 unsigned timeIdx)
const
416 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
417 if (isFineMaterial_(pos))
418 return fineThermalCondParams_;
419 return coarseThermalCondParams_;
437 template <
class Context>
439 const Context& context,
440 unsigned spaceIdx,
unsigned timeIdx)
const
442 const auto& pos = context.cvCenter(spaceIdx, timeIdx);
443 assert(onLeftBoundary_(pos) ||
444 onLowerBoundary_(pos) ||
445 onRightBoundary_(pos) ||
446 onUpperBoundary_(pos));
449 RateVector massRate(0.0);
450 massRate[conti0EqIdx + AirIdx] = -1e-3;
453 values.setMassRate(massRate);
456 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
457 initialFluidState_(fs, context, spaceIdx, timeIdx);
459 Scalar hl = fs.enthalpy(liquidPhaseIdx);
460 Scalar hg = fs.enthalpy(gasPhaseIdx);
461 values.setEnthalpyRate(values[conti0EqIdx + AirIdx] * hg +
462 values[conti0EqIdx + H2OIdx] * hl);
465 else if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
466 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
467 initialFluidState_(fs, context, spaceIdx, timeIdx);
470 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
490 template <
class Context>
492 const Context& context,
494 unsigned timeIdx)
const
496 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
497 initialFluidState_(fs, context, spaceIdx, timeIdx);
500 values.assignMassConservative(fs, matParams,
true);
509 template <
class Context>
511 const Context& context OPM_UNUSED,
512 unsigned spaceIdx OPM_UNUSED,
513 unsigned timeIdx OPM_UNUSED)
const
519 bool onLeftBoundary_(
const GlobalPosition& pos)
const
520 {
return pos[0] < eps_; }
522 bool onRightBoundary_(
const GlobalPosition& pos)
const
523 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
525 bool onLowerBoundary_(
const GlobalPosition& pos)
const
526 {
return pos[1] < eps_; }
528 bool onUpperBoundary_(
const GlobalPosition& pos)
const
529 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
531 bool onInlet_(
const GlobalPosition& pos)
const
532 {
return onLowerBoundary_(pos) && (15.0 < pos[0]) && (pos[0] < 25.0); }
534 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const
535 {
return (20 < pos[0]) && (pos[0] < 30) && (pos[1] < 30); }
537 template <
class Context,
class Flu
idState>
538 void initialFluidState_(FluidState& fs,
539 const Context& context,
541 unsigned timeIdx)
const
543 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
545 Scalar densityW = 1000.0;
546 fs.setPressure(liquidPhaseIdx, 1e5 + (maxDepth_ - pos[1])*densityW*9.81);
547 fs.setSaturation(liquidPhaseIdx, 1.0);
548 fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
549 fs.setMoleFraction(liquidPhaseIdx, AirIdx, 0.0);
551 if (inHighTemperatureRegion_(pos))
552 fs.setTemperature(380);
554 fs.setTemperature(283.0 + (maxDepth_ - pos[1])*0.03);
557 fs.setSaturation(gasPhaseIdx, 0);
558 Scalar pc[numPhases];
560 MaterialLaw::capillaryPressures(pc, matParams, fs);
561 fs.setPressure(gasPhaseIdx, fs.pressure(liquidPhaseIdx) + (pc[gasPhaseIdx] - pc[liquidPhaseIdx]));
563 typename FluidSystem::template ParameterCache<Scalar> paramCache;
564 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
565 CFRP::solve(fs, paramCache, liquidPhaseIdx,
true,
true);
568 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
570 Scalar lambdaGranite = 2.8;
573 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
574 fs.setTemperature(293.15);
575 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
576 fs.setPressure(phaseIdx, 1.0135e5);
579 typename FluidSystem::template ParameterCache<Scalar> paramCache;
580 paramCache.updateAll(fs);
581 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
582 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
583 fs.setDensity(phaseIdx, rho);
586 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
587 Scalar lambdaSaturated;
588 if (FluidSystem::isLiquid(phaseIdx)) {
590 FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
591 lambdaSaturated = std::pow(lambdaGranite, (1-poro)) + std::pow(lambdaFluid, poro);
594 lambdaSaturated = std::pow(lambdaGranite, (1-poro));
596 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
597 if (!FluidSystem::isLiquid(phaseIdx))
598 params.setVacuumLambda(lambdaSaturated);
602 bool isFineMaterial_(
const GlobalPosition& pos)
const
603 {
return pos[dim-1] > layerBottom_; }
609 Scalar finePorosity_;
610 Scalar coarsePorosity_;
612 MaterialLawParams fineMaterialParams_;
613 MaterialLawParams coarseMaterialParams_;
615 ThermalConductionLawParams fineThermalCondParams_;
616 ThermalConductionLawParams coarseThermalCondParams_;
617 SolidEnergyLawParams solidEnergyLawParams_;
Non-isothermal gas injection problem where a air is injected into a fully water saturated medium.
Definition: waterairproblem.hh:206
const SolidEnergyLawParams & solidEnergyLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Return the parameters for the energy storage law of the rock.
Definition: waterairproblem.hh:402
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:359
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:412
void finishInit()
Definition: waterairproblem.hh:268
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:491
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:438
std::string name() const
Definition: waterairproblem.hh:321
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:371
void endTimeStep()
Definition: waterairproblem.hh:334
WaterAirProblem(Simulator &simulator)
Definition: waterairproblem.hh:261
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: waterairproblem.hh:510
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:384
Definition: waterairproblem.hh:64