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>
53 #include <opm/material/common/Unused.hpp>
55 #include <dune/common/version.hh>
56 #include <dune/common/fmatrix.hh>
57 #include <dune/common/fvector.hh>
64 template <
class TypeTag>
65 class FractureProblem;
68 namespace Opm::Properties {
73 struct FractureProblem {
using InheritsFrom = std::tuple<DiscreteFractureModel>; };
77 template<
class TypeTag>
79 {
using type = Dune::ALUGrid<2, 2, Dune::simplex, Dune::nonconforming>; };
82 template<
class TypeTag>
83 struct Vanguard<TypeTag, TTag::
FractureProblem> {
using type = Opm::DgfVanguard<TypeTag>; };
86 template<
class TypeTag>
90 template<
class TypeTag>
94 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
97 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
101 template<
class TypeTag>
105 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
108 using type = Opm::LiquidPhase<Scalar, Opm::DNAPL<Scalar> >;
112 template<
class TypeTag>
116 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
117 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
118 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
120 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
121 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
122 FluidSystem::wettingPhaseIdx,
123 FluidSystem::nonWettingPhaseIdx>;
127 using EffectiveLaw = Opm::RegularizedBrooksCorey<Traits>;
131 using type = Opm::EffToAbsLaw<EffectiveLaw>;
135 template<
class TypeTag>
136 struct EnableEnergy<TypeTag, TTag::
FractureProblem> {
static constexpr
bool value =
true; };
139 template<
class TypeTag>
143 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
144 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
148 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
152 template<
class TypeTag>
154 {
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
157 template<
class TypeTag>
158 struct EnableGravity<TypeTag, TTag::
FractureProblem> {
static constexpr
bool value =
false; };
161 template<
class TypeTag>
162 struct EnableConstraints<TypeTag, TTag::
FractureProblem> {
static constexpr
bool value =
true; };
165 template<
class TypeTag>
166 struct GridFile<TypeTag, TTag::
FractureProblem> {
static constexpr
auto value =
"data/fracture.art.dgf"; };
169 template<
class TypeTag>
172 using type = GetPropType<TypeTag, Scalar>;
173 static constexpr type value = 3e3;
177 template<
class TypeTag>
180 using type = GetPropType<TypeTag, Scalar>;
181 static constexpr type value = 100;
199 template <
class TypeTag>
202 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
203 using GridView = GetPropType<TypeTag, Properties::GridView>;
204 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
205 using WettingPhase = GetPropType<TypeTag, Properties::WettingPhase>;
206 using NonwettingPhase = GetPropType<TypeTag, Properties::NonwettingPhase>;
207 using Constraints = GetPropType<TypeTag, Properties::Constraints>;
208 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
209 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
210 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
211 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
212 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
213 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
214 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
215 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
216 using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
217 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
218 using Model = GetPropType<TypeTag, Properties::Model>;
222 wettingPhaseIdx = MaterialLaw::wettingPhaseIdx,
223 nonWettingPhaseIdx = MaterialLaw::nonWettingPhaseIdx,
226 numPhases = FluidSystem::numPhases,
229 dim = GridView::dimension,
230 dimWorld = GridView::dimensionworld
233 using FluidState = Opm::ImmiscibleFluidState<Scalar, FluidSystem>;
235 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
236 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
241 bool contains(Dune::GeometryType gt)
242 {
return gt.dim() == dim - 1; }
244 using FaceMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
246 using FractureMapper = Opm::FractureMapper<TypeTag>;
253 : ParentType(simulator)
261 ParentType::finishInit();
264 temperature_ = 273.15 + 20;
266 matrixMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.0);
267 matrixMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
268 fractureMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.0);
269 fractureMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
272 matrixMaterialParams_.setEntryPC(0.0);
273 matrixMaterialParams_.setMaxPC(2000.0);
274 fractureMaterialParams_.setEntryPC(0.0);
275 fractureMaterialParams_.setMaxPC(1000.0);
279 matrixMaterialParams_.setEntryPressure(2000);
280 matrixMaterialParams_.setLambda(2.0);
281 matrixMaterialParams_.setPcLowSw(1e-1);
282 fractureMaterialParams_.setEntryPressure(1000);
283 fractureMaterialParams_.setLambda(2.0);
284 fractureMaterialParams_.setPcLowSw(5e-2);
288 matrixMaterialParams_.setVgAlpha(0.0037);
289 matrixMaterialParams_.setVgN(4.7);
290 fractureMaterialParams_.setVgAlpha(0.0025);
291 fractureMaterialParams_.setVgN(4.7);
294 matrixMaterialParams_.finalize();
295 fractureMaterialParams_.finalize();
297 matrixK_ = this->toDimMatrix_(1e-15);
298 fractureK_ = this->toDimMatrix_(1e5 * 1e-15);
300 matrixPorosity_ = 0.10;
301 fracturePorosity_ = 0.25;
302 fractureWidth_ = 1e-3;
305 initEnergyParams_(thermalConductionParams_, matrixPorosity_);
318 std::ostringstream oss;
319 oss <<
"fracture_" << Model::name();
335 this->model().globalStorage(storage);
338 if (this->gridView().comm().rank() == 0) {
339 std::cout <<
"Storage: " << storage << std::endl << std::flush;
347 template <
class Context>
349 unsigned spaceIdx OPM_UNUSED,
350 unsigned timeIdx OPM_UNUSED)
const
351 {
return temperature_; }
363 template <
class Context>
365 unsigned spaceIdx OPM_UNUSED,
366 unsigned timeIdx OPM_UNUSED)
const
374 template <
class Context>
376 unsigned spaceIdx OPM_UNUSED,
377 unsigned timeIdx OPM_UNUSED)
const
378 {
return fractureK_; }
383 template <
class Context>
385 unsigned spaceIdx OPM_UNUSED,
386 unsigned timeIdx OPM_UNUSED)
const
387 {
return matrixPorosity_; }
394 template <
class Context>
396 unsigned spaceIdx OPM_UNUSED,
397 unsigned timeIdx OPM_UNUSED)
const
398 {
return fracturePorosity_; }
403 template <
class Context>
405 unsigned spaceIdx OPM_UNUSED,
406 unsigned timeIdx OPM_UNUSED)
const
407 {
return matrixMaterialParams_; }
414 template <
class Context>
416 unsigned spaceIdx OPM_UNUSED,
417 unsigned timeIdx OPM_UNUSED)
const
418 {
return fractureMaterialParams_; }
424 {
return this->simulator().vanguard().fractureMapper(); }
438 template <
class Context>
440 unsigned spaceIdx1 OPM_UNUSED,
441 unsigned spaceIdx2 OPM_UNUSED,
442 unsigned timeIdx OPM_UNUSED)
const
443 {
return fractureWidth_; }
448 template <
class Context>
449 const ThermalConductionLawParams&
451 unsigned spaceIdx OPM_UNUSED,
452 unsigned timeIdx OPM_UNUSED)
const
453 {
return thermalConductionParams_; }
460 template <
class Context>
461 const SolidEnergyLawParams&
463 unsigned spaceIdx OPM_UNUSED,
464 unsigned timeIdx OPM_UNUSED)
const
465 {
return solidEnergyParams_; }
477 template <
class Context>
478 void boundary(BoundaryRateVector& values,
const Context& context,
479 unsigned spaceIdx,
unsigned timeIdx)
const
481 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
483 if (onRightBoundary_(pos)) {
486 FluidState fluidState;
487 fluidState.setTemperature(temperature_);
489 fluidState.setSaturation(wettingPhaseIdx, 0.0);
490 fluidState.setSaturation(nonWettingPhaseIdx,
491 1.0 - fluidState.saturation(wettingPhaseIdx));
493 fluidState.setPressure(wettingPhaseIdx, 1e5);
494 fluidState.setPressure(nonWettingPhaseIdx, fluidState.pressure(wettingPhaseIdx));
496 typename FluidSystem::template ParameterCache<Scalar> paramCache;
497 paramCache.updateAll(fluidState);
498 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
499 fluidState.setDensity(phaseIdx,
500 FluidSystem::density(fluidState, paramCache, phaseIdx));
501 fluidState.setViscosity(phaseIdx,
502 FluidSystem::viscosity(fluidState, paramCache, phaseIdx));
506 values.setFreeFlow(context, spaceIdx, timeIdx, fluidState);
524 template <
class Context>
526 unsigned spaceIdx,
unsigned timeIdx)
const
528 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
530 if (!onLeftBoundary_(pos))
534 unsigned globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
544 FluidState fractureFluidState;
545 fractureFluidState.setTemperature(temperature_ + 10.0);
547 fractureFluidState.setSaturation(wettingPhaseIdx, 1.0);
548 fractureFluidState.setSaturation(nonWettingPhaseIdx,
549 1.0 - fractureFluidState.saturation(
552 Scalar pCFracture[numPhases];
553 MaterialLaw::capillaryPressures(pCFracture, fractureMaterialParams_,
556 fractureFluidState.setPressure(wettingPhaseIdx, 1.0e5);
557 fractureFluidState.setPressure(nonWettingPhaseIdx,
558 fractureFluidState.pressure(wettingPhaseIdx)
559 + (pCFracture[nonWettingPhaseIdx]
560 - pCFracture[wettingPhaseIdx]));
563 constraints.assignNaiveFromFracture(fractureFluidState,
564 matrixMaterialParams_);
570 template <
class Context>
572 const Context& context OPM_UNUSED,
573 unsigned spaceIdx OPM_UNUSED,
574 unsigned timeIdx OPM_UNUSED)
const
576 FluidState fluidState;
577 fluidState.setTemperature(temperature_);
578 fluidState.setPressure(FluidSystem::wettingPhaseIdx, 1e5);
579 fluidState.setPressure(nonWettingPhaseIdx, fluidState.pressure(wettingPhaseIdx));
581 fluidState.setSaturation(wettingPhaseIdx, 0.0);
582 fluidState.setSaturation(nonWettingPhaseIdx,
583 1.0 - fluidState.saturation(wettingPhaseIdx));
585 values.assignNaive(fluidState);
594 template <
class Context>
596 const Context& context OPM_UNUSED,
597 unsigned spaceIdx OPM_UNUSED,
598 unsigned timeIdx OPM_UNUSED)
const
599 { rate = Scalar(0.0); }
604 bool onLeftBoundary_(
const GlobalPosition& pos)
const
605 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
607 bool onRightBoundary_(
const GlobalPosition& pos)
const
608 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
610 bool onLowerBoundary_(
const GlobalPosition& pos)
const
611 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
613 bool onUpperBoundary_(
const GlobalPosition& pos)
const
614 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
616 void initEnergyParams_(ThermalConductionLawParams& params, Scalar poro)
619 solidEnergyParams_.setSolidHeatCapacity(790.0
621 solidEnergyParams_.finalize();
623 Scalar lambdaGranite = 2.8;
626 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
627 fs.setTemperature(293.15);
628 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
629 fs.setPressure(phaseIdx, 1.0135e5);
632 typename FluidSystem::template ParameterCache<Scalar> paramCache;
633 paramCache.updateAll(fs);
634 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
635 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
636 fs.setDensity(phaseIdx, rho);
639 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
640 Scalar lambdaSaturated;
641 if (FluidSystem::isLiquid(phaseIdx)) {
642 Scalar lambdaFluid = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
644 std::pow(lambdaGranite, (1 - poro))
645 + std::pow(lambdaFluid, poro);
648 lambdaSaturated = std::pow(lambdaGranite, (1 - poro));
650 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
653 Scalar lambdaVac = std::pow(lambdaGranite, (1 - poro));
654 params.setVacuumLambda(lambdaVac);
658 DimMatrix fractureK_;
660 Scalar matrixPorosity_;
661 Scalar fracturePorosity_;
663 Scalar fractureWidth_;
665 MaterialLawParams fractureMaterialParams_;
666 MaterialLawParams matrixMaterialParams_;
668 ThermalConductionLawParams thermalConductionParams_;
669 SolidEnergyLawParams solidEnergyParams_;
Two-phase problem which involves fractures.
Definition: fractureproblem.hh:201
void initial(PrimaryVariables &values, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:571
const MaterialLawParams & materialLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:404
const MaterialLawParams & fractureMaterialLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
The parameters for the material law inside the fractures.
Definition: fractureproblem.hh:415
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:348
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:450
void finishInit()
Definition: fractureproblem.hh:259
std::string name() const
Definition: fractureproblem.hh:316
FractureProblem(Simulator &simulator)
Definition: fractureproblem.hh:252
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:525
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:595
const DimMatrix & intrinsicPermeability(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:364
const DimMatrix & fractureIntrinsicPermeability(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Intrinsic permeability of fractures.
Definition: fractureproblem.hh:375
void endTimeStep()
Called directly after the time integration.
Definition: fractureproblem.hh:326
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: fractureproblem.hh:462
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:384
const FractureMapper & fractureMapper() const
Returns the object representating the fracture topology.
Definition: fractureproblem.hh:423
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:478
Scalar fracturePorosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
The porosity inside the fractures.
Definition: fractureproblem.hh:395
Scalar fractureWidth(const Context &context OPM_UNUSED, unsigned spaceIdx1 OPM_UNUSED, unsigned spaceIdx2 OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Returns the width of the fracture.
Definition: fractureproblem.hh:439
Definition: fractureproblem.hh:73