28#ifndef EWOMS_OBSTACLE_PROBLEM_HH
29#define EWOMS_OBSTACLE_PROBLEM_HH
31#include <dune/common/fmatrix.hh>
32#include <dune/common/fvector.hh>
33#include <dune/common/version.hh>
35#include <dune/grid/yaspgrid.hh>
36#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
38#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
39#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
40#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
41#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
42#include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
43#include <opm/material/fluidstates/CompositionalFluidState.hpp>
44#include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
45#include <opm/material/thermal/ConstantSolidHeatCapLaw.hpp>
46#include <opm/material/thermal/SomertonThermalConductionLaw.hpp>
57template <
class TypeTag>
68template<
class TypeTag>
69struct Grid<TypeTag, TTag::ObstacleBaseProblem> {
using type = Dune::YaspGrid<2>; };
72template<
class TypeTag>
76template<
class TypeTag>
78{
using type = Opm::H2ON2FluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
81template<
class TypeTag>
88 using MaterialTraits = Opm::TwoPhaseMaterialTraits<Scalar,
89 FluidSystem::liquidPhaseIdx,
90 FluidSystem::gasPhaseIdx>;
92 using EffMaterialLaw = Opm::LinearMaterial<MaterialTraits>;
95 using type = Opm::EffToAbsLaw<EffMaterialLaw>;
99template<
class TypeTag>
108 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
112template<
class TypeTag>
114{
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
145template <
class TypeTag>
164 dim = GridView::dimension,
165 dimWorld = GridView::dimensionworld,
166 numPhases = getPropValue<TypeTag, Properties::NumPhases>(),
167 gasPhaseIdx = FluidSystem::gasPhaseIdx,
168 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
169 H2OIdx = FluidSystem::H2OIdx,
170 N2Idx = FluidSystem::N2Idx
173 using GlobalPosition = Dune::FieldVector<typename GridView::ctype, dimWorld>;
174 using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
175 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
184 : ParentType(simulator)
192 ParentType::finishInit();
195 temperature_ = 273.15 + 25;
198 Scalar Tmin = temperature_ - 1.0;
199 Scalar Tmax = temperature_ + 1.0;
202 Scalar pmin = 1.0e5 * 0.75;
203 Scalar pmax = 2.0e5 * 1.25;
206 FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
209 coarseK_ = this->toDimMatrix_(1e-12);
210 fineK_ = this->toDimMatrix_(1e-15);
214 coarsePorosity_ = 0.3;
217 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
218 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
219 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
220 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
224 fineMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
225 fineMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
226 coarseMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
227 coarseMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
239 fineMaterialParams_.finalize();
240 coarseMaterialParams_.finalize();
243 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
244 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
247 solidEnergyLawParams_.setSolidHeatCapacity(790.0
249 solidEnergyLawParams_.finalize();
259 ParentType::registerParameters();
261 Parameters::SetDefault<Parameters::GridFile>(
"./data/obstacle_24x16.dgf");
262 Parameters::SetDefault<Parameters::EndTime<Scalar>>(1e4);
263 Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(250);
264 Parameters::SetDefault<Parameters::EnableGravity>(
true);
273 this->model().checkConservativeness();
276 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
277 PrimaryVariables phaseStorage;
278 this->model().globalPhaseStorage(phaseStorage, phaseIdx);
280 if (this->gridView().comm().rank() == 0) {
281 std::cout <<
"Storage in " << FluidSystem::phaseName(phaseIdx)
282 <<
"Phase: [" << phaseStorage <<
"]"
283 <<
"\n" << std::flush;
289 this->model().globalStorage(storage);
292 if (this->gridView().comm().rank() == 0) {
293 std::cout <<
"Storage total: [" << storage <<
"]"
294 <<
"\n" << std::flush;
309 std::ostringstream oss;
311 <<
"_" << Model::name();
320 template <
class Context>
324 {
return temperature_; }
329 template <
class Context>
333 unsigned timeIdx)
const
335 if (isFineMaterial_(context.pos(spaceIdx, timeIdx)))
343 template <
class Context>
346 unsigned timeIdx)
const
348 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
349 if (isFineMaterial_(pos))
350 return finePorosity_;
352 return coarsePorosity_;
358 template <
class Context>
359 const MaterialLawParams&
362 unsigned timeIdx)
const
364 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
365 if (isFineMaterial_(pos))
366 return fineMaterialParams_;
368 return coarseMaterialParams_;
376 template <
class Context>
377 const SolidEnergyLawParams&
381 {
return solidEnergyLawParams_; }
386 template <
class Context>
387 const ThermalConductionLawParams &
390 unsigned timeIdx)
const
392 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
393 if (isFineMaterial_(pos))
394 return fineThermalCondParams_;
395 return coarseThermalCondParams_;
408 template <
class Context>
410 const Context& context,
412 unsigned timeIdx)
const
414 const auto& pos = context.pos(spaceIdx, timeIdx);
417 values.setFreeFlow(context, spaceIdx, timeIdx, inletFluidState_);
418 else if (onOutlet_(pos))
419 values.setFreeFlow(context, spaceIdx, timeIdx, outletFluidState_);
434 template <
class Context>
436 const Context& context,
438 unsigned timeIdx)
const
441 values.assignMassConservative(outletFluidState_, matParams);
450 template <
class Context>
464 bool isFineMaterial_(
const GlobalPosition& pos)
const
465 {
return 10 <= pos[0] && pos[0] <= 20 && 0 <= pos[1] && pos[1] <= 35; }
467 bool onInlet_(
const GlobalPosition& globalPos)
const
469 Scalar x = globalPos[0];
470 Scalar y = globalPos[1];
471 return x >= 60 - eps_ && y <= 10;
474 bool onOutlet_(
const GlobalPosition& globalPos)
const
476 Scalar x = globalPos[0];
477 Scalar y = globalPos[1];
478 return x < eps_ && y <= 10;
481 void initFluidStates_()
483 initFluidState_(inletFluidState_, coarseMaterialParams_,
485 initFluidState_(outletFluidState_, coarseMaterialParams_,
489 template <
class Flu
idState>
490 void initFluidState_(FluidState& fs,
const MaterialLawParams& matParams,
bool isInlet)
492 unsigned refPhaseIdx;
493 unsigned otherPhaseIdx;
496 fs.setTemperature(temperature_);
500 refPhaseIdx = liquidPhaseIdx;
501 otherPhaseIdx = gasPhaseIdx;
504 fs.setSaturation(liquidPhaseIdx, 1.0);
507 fs.setPressure(liquidPhaseIdx, 2e5);
510 fs.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
511 fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
515 refPhaseIdx = gasPhaseIdx;
516 otherPhaseIdx = liquidPhaseIdx;
519 fs.setSaturation(gasPhaseIdx, 1.0);
522 fs.setPressure(gasPhaseIdx, 1e5);
525 fs.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
526 fs.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
530 fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
534 MaterialLaw::capillaryPressures(pC, matParams, fs);
535 fs.setPressure(otherPhaseIdx, fs.pressure(refPhaseIdx)
536 + (pC[otherPhaseIdx] - pC[refPhaseIdx]));
540 using ComputeFromReferencePhase = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
542 typename FluidSystem::template ParameterCache<Scalar> paramCache;
543 ComputeFromReferencePhase::solve(fs, paramCache, refPhaseIdx,
548 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
550 Scalar lambdaWater = 0.6;
551 Scalar lambdaGranite = 2.8;
553 Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
554 * std::pow(lambdaWater, poro);
555 Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
557 params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
558 params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
559 params.setVacuumLambda(lambdaDry);
565 Scalar coarsePorosity_;
566 Scalar finePorosity_;
568 MaterialLawParams fineMaterialParams_;
569 MaterialLawParams coarseMaterialParams_;
571 ThermalConductionLawParams fineThermalCondParams_;
572 ThermalConductionLawParams coarseThermalCondParams_;
573 SolidEnergyLawParams solidEnergyLawParams_;
575 Opm::CompositionalFluidState<Scalar, FluidSystem> inletFluidState_;
576 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:147
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:344
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: obstacleproblem.hh:190
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:360
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:331
const ThermalConductionLawParams & thermalConductionParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:388
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: obstacleproblem.hh:321
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: obstacleproblem.hh:409
std::string name() const
The problem name.
Definition: obstacleproblem.hh:307
void endTimeStep()
Called by the simulator after each time integration.
Definition: obstacleproblem.hh:270
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: obstacleproblem.hh:435
const SolidEnergyLawParams & solidEnergyLawParams(const Context &, unsigned, unsigned) const
Return the parameters for the energy storage law of the rock.
Definition: obstacleproblem.hh:378
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: obstacleproblem.hh:451
ObstacleProblem(Simulator &simulator)
Definition: obstacleproblem.hh:183
static void registerParameters()
Definition: obstacleproblem.hh:257
Defines the common parameters for the porous medium multi-phase models.
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:235
Declares the properties required for the NCP compositional multi-phase model.
Opm::H2ON2FluidSystem< GetPropType< TypeTag, Properties::Scalar > > type
Definition: obstacleproblem.hh:78
The fluid systems including the information about the phases.
Definition: multiphasebaseproperties.hh:69
Dune::YaspGrid< 2 > type
Definition: obstacleproblem.hh:69
The type of the DUNE grid.
Definition: basicproperties.hh:100
Opm::EffToAbsLaw< EffMaterialLaw > type
Definition: obstacleproblem.hh:95
The material law which ought to be used (extracted from the spatial parameters)
Definition: multiphasebaseproperties.hh:51
The type of the problem.
Definition: fvbaseproperties.hh:81
Opm::ConstantSolidHeatCapLaw< GetPropType< TypeTag, Properties::Scalar > > type
Definition: obstacleproblem.hh:114
The material law for the energy stored in the solid matrix.
Definition: multiphasebaseproperties.hh:57
Definition: obstacleproblem.hh:64
Opm::SomertonThermalConductionLaw< FluidSystem, Scalar > type
Definition: obstacleproblem.hh:108
The material law for thermal conduction.
Definition: multiphasebaseproperties.hh:63