28#ifndef EWOMS_FRACTURE_PROBLEM_HH
29#define EWOMS_FRACTURE_PROBLEM_HH
33#define DISABLE_ALUGRID_SFC_ORDERING 1
34#include <dune/alugrid/grid.hh>
35#include <dune/alugrid/dgf.hh>
37#error "dune-alugrid not found!"
40#include <opm/models/discretefracture/discretefracturemodel.hh>
41#include <opm/models/io/dgfvanguard.hh>
43#include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
44#include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
45#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
46#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
47#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
48#include <opm/material/thermal/SomertonThermalConductionLaw.hpp>
49#include <opm/material/thermal/ConstantSolidHeatCapLaw.hpp>
50#include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp>
51#include <opm/material/components/SimpleH2O.hpp>
52#include <opm/material/components/Dnapl.hpp>
54#include <dune/common/version.hh>
55#include <dune/common/fmatrix.hh>
56#include <dune/common/fvector.hh>
63template <
class TypeTag>
67namespace Opm::Properties {
72struct FractureProblem {
using InheritsFrom = std::tuple<DiscreteFractureModel>; };
76template<
class TypeTag>
78{
using type = Dune::ALUGrid<2, 2, Dune::simplex, Dune::nonconforming>; };
81template<
class TypeTag>
82struct Vanguard<TypeTag, TTag::
FractureProblem> {
using type = Opm::DgfVanguard<TypeTag>; };
85template<
class TypeTag>
89template<
class TypeTag>
93 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
96 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
100template<
class TypeTag>
104 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
107 using type = Opm::LiquidPhase<Scalar, Opm::DNAPL<Scalar> >;
111template<
class TypeTag>
115 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
116 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
117 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
119 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
120 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
121 FluidSystem::wettingPhaseIdx,
122 FluidSystem::nonWettingPhaseIdx>;
126 using EffectiveLaw = Opm::RegularizedBrooksCorey<Traits>;
130 using type = Opm::EffToAbsLaw<EffectiveLaw>;
134template<
class TypeTag>
135struct EnableEnergy<TypeTag, TTag::
FractureProblem> {
static constexpr bool value =
true; };
138template<
class TypeTag>
142 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
143 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
147 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
151template<
class TypeTag>
153{
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
156template<
class TypeTag>
157struct EnableGravity<TypeTag, TTag::
FractureProblem> {
static constexpr bool value =
false; };
160template<
class TypeTag>
161struct EnableConstraints<TypeTag, TTag::
FractureProblem> {
static constexpr bool value =
true; };
164template<
class TypeTag>
165struct GridFile<TypeTag, TTag::
FractureProblem> {
static constexpr auto value =
"data/fracture.art.dgf"; };
168template<
class TypeTag>
171 using type = GetPropType<TypeTag, Scalar>;
172 static constexpr type value = 3e3;
176template<
class TypeTag>
179 using type = GetPropType<TypeTag, Scalar>;
180 static constexpr type value = 100;
198template <
class TypeTag>
201 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
202 using GridView = GetPropType<TypeTag, Properties::GridView>;
203 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
204 using WettingPhase = GetPropType<TypeTag, Properties::WettingPhase>;
205 using NonwettingPhase = GetPropType<TypeTag, Properties::NonwettingPhase>;
206 using Constraints = GetPropType<TypeTag, Properties::Constraints>;
207 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
208 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
209 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
210 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
211 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
212 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
213 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
214 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
215 using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
216 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
217 using Model = GetPropType<TypeTag, Properties::Model>;
221 wettingPhaseIdx = MaterialLaw::wettingPhaseIdx,
222 nonWettingPhaseIdx = MaterialLaw::nonWettingPhaseIdx,
225 numPhases = FluidSystem::numPhases,
228 dim = GridView::dimension,
229 dimWorld = GridView::dimensionworld
232 using FluidState = Opm::ImmiscibleFluidState<Scalar, FluidSystem>;
234 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
235 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
240 bool contains(Dune::GeometryType gt)
241 {
return gt.dim() == dim - 1; }
243 using FaceMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
245 using FractureMapper = Opm::FractureMapper<TypeTag>;
252 : ParentType(simulator)
260 ParentType::finishInit();
263 temperature_ = 273.15 + 20;
265 matrixMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.0);
266 matrixMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
267 fractureMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.0);
268 fractureMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
271 matrixMaterialParams_.setEntryPC(0.0);
272 matrixMaterialParams_.setMaxPC(2000.0);
273 fractureMaterialParams_.setEntryPC(0.0);
274 fractureMaterialParams_.setMaxPC(1000.0);
278 matrixMaterialParams_.setEntryPressure(2000);
279 matrixMaterialParams_.setLambda(2.0);
280 matrixMaterialParams_.setPcLowSw(1e-1);
281 fractureMaterialParams_.setEntryPressure(1000);
282 fractureMaterialParams_.setLambda(2.0);
283 fractureMaterialParams_.setPcLowSw(5e-2);
287 matrixMaterialParams_.setVgAlpha(0.0037);
288 matrixMaterialParams_.setVgN(4.7);
289 fractureMaterialParams_.setVgAlpha(0.0025);
290 fractureMaterialParams_.setVgN(4.7);
293 matrixMaterialParams_.finalize();
294 fractureMaterialParams_.finalize();
296 matrixK_ = this->toDimMatrix_(1e-15);
297 fractureK_ = this->toDimMatrix_(1e5 * 1e-15);
299 matrixPorosity_ = 0.10;
300 fracturePorosity_ = 0.25;
301 fractureWidth_ = 1e-3;
304 initEnergyParams_(thermalConductionParams_, matrixPorosity_);
317 std::ostringstream oss;
318 oss <<
"fracture_" << Model::name();
334 this->model().globalStorage(storage);
337 if (this->gridView().comm().rank() == 0) {
338 std::cout <<
"Storage: " << storage << std::endl << std::flush;
346 template <
class Context>
348 [[maybe_unused]]
unsigned spaceIdx,
349 [[maybe_unused]]
unsigned timeIdx)
const
350 {
return temperature_; }
362 template <
class Context>
364 [[maybe_unused]]
unsigned spaceIdx,
365 [[maybe_unused]]
unsigned timeIdx)
const
373 template <
class Context>
375 [[maybe_unused]]
unsigned spaceIdx,
376 [[maybe_unused]]
unsigned timeIdx)
const
377 {
return fractureK_; }
382 template <
class Context>
383 Scalar
porosity([[maybe_unused]]
const Context& context,
384 [[maybe_unused]]
unsigned spaceIdx,
385 [[maybe_unused]]
unsigned timeIdx)
const
386 {
return matrixPorosity_; }
393 template <
class Context>
395 [[maybe_unused]]
unsigned spaceIdx,
396 [[maybe_unused]]
unsigned timeIdx)
const
397 {
return fracturePorosity_; }
402 template <
class Context>
404 [[maybe_unused]]
unsigned spaceIdx,
405 [[maybe_unused]]
unsigned timeIdx)
const
406 {
return matrixMaterialParams_; }
413 template <
class Context>
415 [[maybe_unused]]
unsigned spaceIdx,
416 [[maybe_unused]]
unsigned timeIdx)
const
417 {
return fractureMaterialParams_; }
423 {
return this->simulator().vanguard().fractureMapper(); }
437 template <
class Context>
439 [[maybe_unused]]
unsigned spaceIdx1,
440 [[maybe_unused]]
unsigned spaceIdx2,
441 [[maybe_unused]]
unsigned timeIdx)
const
442 {
return fractureWidth_; }
447 template <
class Context>
448 const ThermalConductionLawParams&
450 [[maybe_unused]]
unsigned spaceIdx,
451 [[maybe_unused]]
unsigned timeIdx)
const
452 {
return thermalConductionParams_; }
459 template <
class Context>
460 const SolidEnergyLawParams&
462 [[maybe_unused]]
unsigned spaceIdx,
463 [[maybe_unused]]
unsigned timeIdx)
const
464 {
return solidEnergyParams_; }
476 template <
class Context>
477 void boundary(BoundaryRateVector& values,
const Context& context,
478 unsigned spaceIdx,
unsigned timeIdx)
const
480 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
482 if (onRightBoundary_(pos)) {
485 FluidState fluidState;
486 fluidState.setTemperature(temperature_);
488 fluidState.setSaturation(wettingPhaseIdx, 0.0);
489 fluidState.setSaturation(nonWettingPhaseIdx,
490 1.0 - fluidState.saturation(wettingPhaseIdx));
492 fluidState.setPressure(wettingPhaseIdx, 1e5);
493 fluidState.setPressure(nonWettingPhaseIdx, fluidState.pressure(wettingPhaseIdx));
495 typename FluidSystem::template ParameterCache<Scalar> paramCache;
496 paramCache.updateAll(fluidState);
497 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
498 fluidState.setDensity(phaseIdx,
499 FluidSystem::density(fluidState, paramCache, phaseIdx));
500 fluidState.setViscosity(phaseIdx,
501 FluidSystem::viscosity(fluidState, paramCache, phaseIdx));
505 values.setFreeFlow(context, spaceIdx, timeIdx, fluidState);
523 template <
class Context>
525 unsigned spaceIdx,
unsigned timeIdx)
const
527 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
529 if (!onLeftBoundary_(pos))
533 unsigned globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
543 FluidState fractureFluidState;
544 fractureFluidState.setTemperature(temperature_ + 10.0);
546 fractureFluidState.setSaturation(wettingPhaseIdx, 1.0);
547 fractureFluidState.setSaturation(nonWettingPhaseIdx,
548 1.0 - fractureFluidState.saturation(
551 Scalar pCFracture[numPhases];
552 MaterialLaw::capillaryPressures(pCFracture, fractureMaterialParams_,
555 fractureFluidState.setPressure(wettingPhaseIdx, 1.0e5);
556 fractureFluidState.setPressure(nonWettingPhaseIdx,
557 fractureFluidState.pressure(wettingPhaseIdx)
558 + (pCFracture[nonWettingPhaseIdx]
559 - pCFracture[wettingPhaseIdx]));
562 constraints.assignNaiveFromFracture(fractureFluidState,
563 matrixMaterialParams_);
569 template <
class Context>
571 [[maybe_unused]]
const Context& context,
572 [[maybe_unused]]
unsigned spaceIdx,
573 [[maybe_unused]]
unsigned timeIdx)
const
575 FluidState fluidState;
576 fluidState.setTemperature(temperature_);
577 fluidState.setPressure(FluidSystem::wettingPhaseIdx, 1e5);
578 fluidState.setPressure(nonWettingPhaseIdx, fluidState.pressure(wettingPhaseIdx));
580 fluidState.setSaturation(wettingPhaseIdx, 0.0);
581 fluidState.setSaturation(nonWettingPhaseIdx,
582 1.0 - fluidState.saturation(wettingPhaseIdx));
584 values.assignNaive(fluidState);
593 template <
class Context>
595 [[maybe_unused]]
const Context& context,
596 [[maybe_unused]]
unsigned spaceIdx,
597 [[maybe_unused]]
unsigned timeIdx)
const
598 { rate = Scalar(0.0); }
603 bool onLeftBoundary_(
const GlobalPosition& pos)
const
604 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
606 bool onRightBoundary_(
const GlobalPosition& pos)
const
607 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
609 bool onLowerBoundary_(
const GlobalPosition& pos)
const
610 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
612 bool onUpperBoundary_(
const GlobalPosition& pos)
const
613 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
615 void initEnergyParams_(ThermalConductionLawParams& params, Scalar poro)
618 solidEnergyParams_.setSolidHeatCapacity(790.0
620 solidEnergyParams_.finalize();
622 Scalar lambdaGranite = 2.8;
625 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
626 fs.setTemperature(293.15);
627 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
628 fs.setPressure(phaseIdx, 1.0135e5);
631 typename FluidSystem::template ParameterCache<Scalar> paramCache;
632 paramCache.updateAll(fs);
633 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
634 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
635 fs.setDensity(phaseIdx, rho);
638 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
639 Scalar lambdaSaturated;
640 if (FluidSystem::isLiquid(phaseIdx)) {
641 Scalar lambdaFluid = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
643 std::pow(lambdaGranite, (1 - poro))
644 + std::pow(lambdaFluid, poro);
647 lambdaSaturated = std::pow(lambdaGranite, (1 - poro));
649 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
652 Scalar lambdaVac = std::pow(lambdaGranite, (1 - poro));
653 params.setVacuumLambda(lambdaVac);
657 DimMatrix fractureK_;
659 Scalar matrixPorosity_;
660 Scalar fracturePorosity_;
662 Scalar fractureWidth_;
664 MaterialLawParams fractureMaterialParams_;
665 MaterialLawParams matrixMaterialParams_;
667 ThermalConductionLawParams thermalConductionParams_;
668 SolidEnergyLawParams solidEnergyParams_;
Two-phase problem which involves fractures.
Definition: fractureproblem.hh:200
const SolidEnergyLawParams & solidEnergyLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the parameters for the energy storage law of the rock.
Definition: fractureproblem.hh:461
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:363
Scalar fractureWidth(const Context &context, unsigned spaceIdx1, unsigned spaceIdx2, unsigned timeIdx) const
Returns the width of the fracture.
Definition: fractureproblem.hh:438
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:570
Scalar fracturePorosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
The porosity inside the fractures.
Definition: fractureproblem.hh:394
void finishInit()
Definition: fractureproblem.hh:258
std::string name() const
Definition: fractureproblem.hh:315
FractureProblem(Simulator &simulator)
Definition: fractureproblem.hh:251
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:524
const DimMatrix & fractureIntrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Intrinsic permeability of fractures.
Definition: fractureproblem.hh:374
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:403
void endTimeStep()
Called directly after the time integration.
Definition: fractureproblem.hh:325
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:449
void source(RateVector &rate, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:594
const MaterialLawParams & fractureMaterialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
The parameters for the material law inside the fractures.
Definition: fractureproblem.hh:414
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:347
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:477
const FractureMapper & fractureMapper() const
Returns the object representating the fracture topology.
Definition: fractureproblem.hh:422
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:383
Definition: fractureproblem.hh:72