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>
42 #include <opm/material/common/Unused.hpp>
44 #include <dune/grid/yaspgrid.hh>
45 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
47 #include <dune/common/version.hh>
48 #include <dune/common/fvector.hh>
49 #include <dune/common/fmatrix.hh>
56 template <
class TypeTag>
57 class ObstacleProblem;
60 namespace Opm::Properties {
67 template<
class TypeTag>
68 struct Grid<TypeTag, TTag::ObstacleBaseProblem> {
using type = Dune::YaspGrid<2>; };
71 template<
class TypeTag>
75 template<
class TypeTag>
76 struct FluidSystem<TypeTag, TTag::ObstacleBaseProblem>
77 {
using type = Opm::H2ON2FluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
80 template<
class TypeTag>
81 struct MaterialLaw<TypeTag, TTag::ObstacleBaseProblem>
85 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
86 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
87 using MaterialTraits = Opm::TwoPhaseMaterialTraits<Scalar,
88 FluidSystem::liquidPhaseIdx,
89 FluidSystem::gasPhaseIdx>;
91 using EffMaterialLaw = Opm::LinearMaterial<MaterialTraits>;
94 using type = Opm::EffToAbsLaw<EffMaterialLaw>;
98 template<
class TypeTag>
99 struct ThermalConductionLaw<TypeTag, TTag::ObstacleBaseProblem>
102 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
103 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
107 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
111 template<
class TypeTag>
112 struct SolidEnergyLaw<TypeTag, TTag::ObstacleBaseProblem>
113 {
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
116 template<
class TypeTag>
117 struct EnableGravity<TypeTag, TTag::ObstacleBaseProblem> {
static constexpr
bool value =
true; };
120 template<
class TypeTag>
121 struct EndTime<TypeTag, TTag::ObstacleBaseProblem>
123 using type = GetPropType<TypeTag, Scalar>;
124 static constexpr type value = 1e4;
128 template<
class TypeTag>
129 struct InitialTimeStepSize<TypeTag, TTag::ObstacleBaseProblem>
131 using type = GetPropType<TypeTag, Scalar>;
132 static constexpr type value = 250;
136 template<
class TypeTag>
137 struct GridFile<TypeTag, TTag::ObstacleBaseProblem> {
static constexpr
auto value =
"./data/obstacle_24x16.dgf"; };
168 template <
class TypeTag>
171 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
173 using GridView = GetPropType<TypeTag, Properties::GridView>;
174 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
175 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
176 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
177 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
178 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
179 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
180 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
181 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
182 using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
183 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
187 dim = GridView::dimension,
188 dimWorld = GridView::dimensionworld,
189 numPhases = getPropValue<TypeTag, Properties::NumPhases>(),
190 gasPhaseIdx = FluidSystem::gasPhaseIdx,
191 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
192 H2OIdx = FluidSystem::H2OIdx,
193 N2Idx = FluidSystem::N2Idx
196 using GlobalPosition = Dune::FieldVector<typename GridView::ctype, dimWorld>;
197 using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
198 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
199 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
200 using Model = GetPropType<TypeTag, Properties::Model>;
207 : ParentType(simulator)
215 ParentType::finishInit();
218 temperature_ = 273.15 + 25;
221 Scalar Tmin = temperature_ - 1.0;
222 Scalar Tmax = temperature_ + 1.0;
225 Scalar pmin = 1.0e5 * 0.75;
226 Scalar pmax = 2.0e5 * 1.25;
229 FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
232 coarseK_ = this->toDimMatrix_(1e-12);
233 fineK_ = this->toDimMatrix_(1e-15);
237 coarsePorosity_ = 0.3;
240 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
241 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
242 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
243 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
247 fineMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
248 fineMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
249 coarseMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
250 coarseMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
262 fineMaterialParams_.finalize();
263 coarseMaterialParams_.finalize();
266 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
267 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
270 solidEnergyLawParams_.setSolidHeatCapacity(790.0
272 solidEnergyLawParams_.finalize();
283 this->model().checkConservativeness();
286 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
287 PrimaryVariables phaseStorage;
288 this->model().globalPhaseStorage(phaseStorage, phaseIdx);
290 if (this->gridView().comm().rank() == 0) {
291 std::cout <<
"Storage in " << FluidSystem::phaseName(phaseIdx)
292 <<
"Phase: [" << phaseStorage <<
"]"
293 <<
"\n" << std::flush;
299 this->model().globalStorage(storage);
302 if (this->gridView().comm().rank() == 0) {
303 std::cout <<
"Storage total: [" << storage <<
"]"
304 <<
"\n" << std::flush;
319 std::ostringstream oss;
321 <<
"_" << Model::name();
330 template <
class Context>
332 unsigned spaceIdx OPM_UNUSED,
333 unsigned timeIdx OPM_UNUSED)
const
334 {
return temperature_; }
339 template <
class Context>
343 unsigned timeIdx)
const
345 if (isFineMaterial_(context.pos(spaceIdx, timeIdx)))
353 template <
class Context>
356 unsigned timeIdx)
const
358 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
359 if (isFineMaterial_(pos))
360 return finePorosity_;
362 return coarsePorosity_;
368 template <
class Context>
369 const MaterialLawParams&
372 unsigned timeIdx)
const
374 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
375 if (isFineMaterial_(pos))
376 return fineMaterialParams_;
378 return coarseMaterialParams_;
386 template <
class Context>
387 const SolidEnergyLawParams&
389 unsigned spaceIdx OPM_UNUSED,
390 unsigned timeIdx OPM_UNUSED)
const
391 {
return solidEnergyLawParams_; }
396 template <
class Context>
397 const ThermalConductionLawParams &
400 unsigned timeIdx)
const
402 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
403 if (isFineMaterial_(pos))
404 return fineThermalCondParams_;
405 return coarseThermalCondParams_;
418 template <
class Context>
420 const Context& context,
422 unsigned timeIdx)
const
424 const auto& pos = context.pos(spaceIdx, timeIdx);
427 values.setFreeFlow(context, spaceIdx, timeIdx, inletFluidState_);
428 else if (onOutlet_(pos))
429 values.setFreeFlow(context, spaceIdx, timeIdx, outletFluidState_);
444 template <
class Context>
446 const Context& context,
448 unsigned timeIdx)
const
451 values.assignMassConservative(outletFluidState_, matParams);
460 template <
class Context>
462 const Context& context OPM_UNUSED,
463 unsigned spaceIdx OPM_UNUSED,
464 unsigned timeIdx OPM_UNUSED)
const
474 bool isFineMaterial_(
const GlobalPosition& pos)
const
475 {
return 10 <= pos[0] && pos[0] <= 20 && 0 <= pos[1] && pos[1] <= 35; }
477 bool onInlet_(
const GlobalPosition& globalPos)
const
479 Scalar x = globalPos[0];
480 Scalar y = globalPos[1];
481 return x >= 60 - eps_ && y <= 10;
484 bool onOutlet_(
const GlobalPosition& globalPos)
const
486 Scalar x = globalPos[0];
487 Scalar y = globalPos[1];
488 return x < eps_ && y <= 10;
491 void initFluidStates_()
493 initFluidState_(inletFluidState_, coarseMaterialParams_,
495 initFluidState_(outletFluidState_, coarseMaterialParams_,
499 template <
class Flu
idState>
500 void initFluidState_(FluidState& fs,
const MaterialLawParams& matParams,
bool isInlet)
502 unsigned refPhaseIdx;
503 unsigned otherPhaseIdx;
506 fs.setTemperature(temperature_);
510 refPhaseIdx = liquidPhaseIdx;
511 otherPhaseIdx = gasPhaseIdx;
514 fs.setSaturation(liquidPhaseIdx, 1.0);
517 fs.setPressure(liquidPhaseIdx, 2e5);
520 fs.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
521 fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
525 refPhaseIdx = gasPhaseIdx;
526 otherPhaseIdx = liquidPhaseIdx;
529 fs.setSaturation(gasPhaseIdx, 1.0);
532 fs.setPressure(gasPhaseIdx, 1e5);
535 fs.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
536 fs.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
540 fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
544 MaterialLaw::capillaryPressures(pC, matParams, fs);
545 fs.setPressure(otherPhaseIdx, fs.pressure(refPhaseIdx)
546 + (pC[otherPhaseIdx] - pC[refPhaseIdx]));
550 using ComputeFromReferencePhase = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
552 typename FluidSystem::template ParameterCache<Scalar> paramCache;
553 ComputeFromReferencePhase::solve(fs, paramCache, refPhaseIdx,
558 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
560 Scalar lambdaWater = 0.6;
561 Scalar lambdaGranite = 2.8;
563 Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
564 * std::pow(lambdaWater, poro);
565 Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
567 params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
568 params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
569 params.setVacuumLambda(lambdaDry);
575 Scalar coarsePorosity_;
576 Scalar finePorosity_;
578 MaterialLawParams fineMaterialParams_;
579 MaterialLawParams coarseMaterialParams_;
581 ThermalConductionLawParams fineThermalCondParams_;
582 ThermalConductionLawParams coarseThermalCondParams_;
583 SolidEnergyLawParams solidEnergyLawParams_;
585 Opm::CompositionalFluidState<Scalar, FluidSystem> inletFluidState_;
586 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:170
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:354
void finishInit()
Definition: obstacleproblem.hh:213
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:370
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: obstacleproblem.hh:461
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: obstacleproblem.hh:331
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:341
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:419
std::string name() const
Definition: obstacleproblem.hh:317
void endTimeStep()
Definition: obstacleproblem.hh:280
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:445
ObstacleProblem(Simulator &simulator)
Definition: obstacleproblem.hh:206
const ThermalConductionLawParams & thermalConductionParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:398
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: obstacleproblem.hh:388
Definition: obstacleproblem.hh:63