28 #ifndef EWOMS_OBSTACLE_PROBLEM_HH
29 #define EWOMS_OBSTACLE_PROBLEM_HH
31 #include <opm/models/ncp/ncpproperties.hh>
33 #include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
34 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
35 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
36 #include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
37 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
38 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
39 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
40 #include <opm/material/thermal/ConstantSolidHeatCapLaw.hpp>
41 #include <opm/material/thermal/SomertonThermalConductionLaw.hpp>
43 #include <dune/grid/yaspgrid.hh>
44 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
46 #include <dune/common/version.hh>
47 #include <dune/common/fvector.hh>
48 #include <dune/common/fmatrix.hh>
55 template <
class TypeTag>
56 class ObstacleProblem;
59 namespace Opm::Properties {
66 template<
class TypeTag>
67 struct Grid<TypeTag, TTag::ObstacleBaseProblem> {
using type = Dune::YaspGrid<2>; };
70 template<
class TypeTag>
74 template<
class TypeTag>
75 struct FluidSystem<TypeTag, TTag::ObstacleBaseProblem>
76 {
using type = Opm::H2ON2FluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
79 template<
class TypeTag>
80 struct MaterialLaw<TypeTag, TTag::ObstacleBaseProblem>
84 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
85 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
86 using MaterialTraits = Opm::TwoPhaseMaterialTraits<Scalar,
87 FluidSystem::liquidPhaseIdx,
88 FluidSystem::gasPhaseIdx>;
90 using EffMaterialLaw = Opm::LinearMaterial<MaterialTraits>;
93 using type = Opm::EffToAbsLaw<EffMaterialLaw>;
97 template<
class TypeTag>
98 struct ThermalConductionLaw<TypeTag, TTag::ObstacleBaseProblem>
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::ObstacleBaseProblem>
112 {
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
115 template<
class TypeTag>
116 struct EnableGravity<TypeTag, TTag::ObstacleBaseProblem> {
static constexpr
bool value =
true; };
119 template<
class TypeTag>
120 struct EndTime<TypeTag, TTag::ObstacleBaseProblem>
122 using type = GetPropType<TypeTag, Scalar>;
123 static constexpr type value = 1e4;
127 template<
class TypeTag>
128 struct InitialTimeStepSize<TypeTag, TTag::ObstacleBaseProblem>
130 using type = GetPropType<TypeTag, Scalar>;
131 static constexpr type value = 250;
135 template<
class TypeTag>
136 struct GridFile<TypeTag, TTag::ObstacleBaseProblem> {
static constexpr
auto value =
"./data/obstacle_24x16.dgf"; };
167 template <
class TypeTag>
170 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
172 using GridView = GetPropType<TypeTag, Properties::GridView>;
173 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
174 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
175 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
176 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
177 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
178 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
179 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
180 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
181 using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
182 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
186 dim = GridView::dimension,
187 dimWorld = GridView::dimensionworld,
188 numPhases = getPropValue<TypeTag, Properties::NumPhases>(),
189 gasPhaseIdx = FluidSystem::gasPhaseIdx,
190 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
191 H2OIdx = FluidSystem::H2OIdx,
192 N2Idx = FluidSystem::N2Idx
195 using GlobalPosition = Dune::FieldVector<typename GridView::ctype, dimWorld>;
196 using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
197 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
198 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
199 using Model = GetPropType<TypeTag, Properties::Model>;
206 : ParentType(simulator)
214 ParentType::finishInit();
217 temperature_ = 273.15 + 25;
220 Scalar Tmin = temperature_ - 1.0;
221 Scalar Tmax = temperature_ + 1.0;
224 Scalar pmin = 1.0e5 * 0.75;
225 Scalar pmax = 2.0e5 * 1.25;
228 FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
231 coarseK_ = this->toDimMatrix_(1e-12);
232 fineK_ = this->toDimMatrix_(1e-15);
236 coarsePorosity_ = 0.3;
239 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
240 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
241 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
242 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
246 fineMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
247 fineMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
248 coarseMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
249 coarseMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
261 fineMaterialParams_.finalize();
262 coarseMaterialParams_.finalize();
265 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
266 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
269 solidEnergyLawParams_.setSolidHeatCapacity(790.0
271 solidEnergyLawParams_.finalize();
282 this->model().checkConservativeness();
285 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
286 PrimaryVariables phaseStorage;
287 this->model().globalPhaseStorage(phaseStorage, phaseIdx);
289 if (this->gridView().comm().rank() == 0) {
290 std::cout <<
"Storage in " << FluidSystem::phaseName(phaseIdx)
291 <<
"Phase: [" << phaseStorage <<
"]"
292 <<
"\n" << std::flush;
298 this->model().globalStorage(storage);
301 if (this->gridView().comm().rank() == 0) {
302 std::cout <<
"Storage total: [" << storage <<
"]"
303 <<
"\n" << std::flush;
318 std::ostringstream oss;
320 <<
"_" << Model::name();
329 template <
class Context>
333 {
return temperature_; }
338 template <
class Context>
342 unsigned timeIdx)
const
344 if (isFineMaterial_(context.pos(spaceIdx, timeIdx)))
352 template <
class Context>
355 unsigned timeIdx)
const
357 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
358 if (isFineMaterial_(pos))
359 return finePorosity_;
361 return coarsePorosity_;
367 template <
class Context>
368 const MaterialLawParams&
371 unsigned timeIdx)
const
373 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
374 if (isFineMaterial_(pos))
375 return fineMaterialParams_;
377 return coarseMaterialParams_;
385 template <
class Context>
386 const SolidEnergyLawParams&
390 {
return solidEnergyLawParams_; }
395 template <
class Context>
396 const ThermalConductionLawParams &
399 unsigned timeIdx)
const
401 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
402 if (isFineMaterial_(pos))
403 return fineThermalCondParams_;
404 return coarseThermalCondParams_;
417 template <
class Context>
419 const Context& context,
421 unsigned timeIdx)
const
423 const auto& pos = context.pos(spaceIdx, timeIdx);
426 values.setFreeFlow(context, spaceIdx, timeIdx, inletFluidState_);
427 else if (onOutlet_(pos))
428 values.setFreeFlow(context, spaceIdx, timeIdx, outletFluidState_);
443 template <
class Context>
445 const Context& context,
447 unsigned timeIdx)
const
450 values.assignMassConservative(outletFluidState_, matParams);
459 template <
class Context>
473 bool isFineMaterial_(
const GlobalPosition& pos)
const
474 {
return 10 <= pos[0] && pos[0] <= 20 && 0 <= pos[1] && pos[1] <= 35; }
476 bool onInlet_(
const GlobalPosition& globalPos)
const
478 Scalar x = globalPos[0];
479 Scalar y = globalPos[1];
480 return x >= 60 - eps_ && y <= 10;
483 bool onOutlet_(
const GlobalPosition& globalPos)
const
485 Scalar x = globalPos[0];
486 Scalar y = globalPos[1];
487 return x < eps_ && y <= 10;
490 void initFluidStates_()
492 initFluidState_(inletFluidState_, coarseMaterialParams_,
494 initFluidState_(outletFluidState_, coarseMaterialParams_,
498 template <
class Flu
idState>
499 void initFluidState_(FluidState& fs,
const MaterialLawParams& matParams,
bool isInlet)
501 unsigned refPhaseIdx;
502 unsigned otherPhaseIdx;
505 fs.setTemperature(temperature_);
509 refPhaseIdx = liquidPhaseIdx;
510 otherPhaseIdx = gasPhaseIdx;
513 fs.setSaturation(liquidPhaseIdx, 1.0);
516 fs.setPressure(liquidPhaseIdx, 2e5);
519 fs.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
520 fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
524 refPhaseIdx = gasPhaseIdx;
525 otherPhaseIdx = liquidPhaseIdx;
528 fs.setSaturation(gasPhaseIdx, 1.0);
531 fs.setPressure(gasPhaseIdx, 1e5);
534 fs.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
535 fs.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
539 fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
543 MaterialLaw::capillaryPressures(pC, matParams, fs);
544 fs.setPressure(otherPhaseIdx, fs.pressure(refPhaseIdx)
545 + (pC[otherPhaseIdx] - pC[refPhaseIdx]));
549 using ComputeFromReferencePhase = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
551 typename FluidSystem::template ParameterCache<Scalar> paramCache;
552 ComputeFromReferencePhase::solve(fs, paramCache, refPhaseIdx,
557 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
559 Scalar lambdaWater = 0.6;
560 Scalar lambdaGranite = 2.8;
562 Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
563 * std::pow(lambdaWater, poro);
564 Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
566 params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
567 params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
568 params.setVacuumLambda(lambdaDry);
574 Scalar coarsePorosity_;
575 Scalar finePorosity_;
577 MaterialLawParams fineMaterialParams_;
578 MaterialLawParams coarseMaterialParams_;
580 ThermalConductionLawParams fineThermalCondParams_;
581 ThermalConductionLawParams coarseThermalCondParams_;
582 SolidEnergyLawParams solidEnergyLawParams_;
584 Opm::CompositionalFluidState<Scalar, FluidSystem> inletFluidState_;
585 Opm::CompositionalFluidState<Scalar, FluidSystem> outletFluidState_;
Problem where liquid water is first stopped by a low-permeability lens and then seeps though it.
Definition: obstacleproblem.hh:169
const SolidEnergyLawParams & solidEnergyLawParams(const Context &, unsigned, unsigned) const
Return the parameters for the energy storage law of the rock.
Definition: obstacleproblem.hh:387
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:353
void finishInit()
Definition: obstacleproblem.hh:212
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:369
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: obstacleproblem.hh:330
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:340
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:418
std::string name() const
Definition: obstacleproblem.hh:316
void endTimeStep()
Definition: obstacleproblem.hh:279
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:444
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: obstacleproblem.hh:460
ObstacleProblem(Simulator &simulator)
Definition: obstacleproblem.hh:205
const ThermalConductionLawParams & thermalConductionParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:397
Definition: obstacleproblem.hh:62