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!"
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>
76template<
class TypeTag>
81template<
class TypeTag>
85template<
class TypeTag>
89template<
class TypeTag>
96 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
100template<
class TypeTag>
107 using type = Opm::LiquidPhase<Scalar, Opm::DNAPL<Scalar> >;
111template<
class TypeTag>
116 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
117 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
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>
138template<
class TypeTag>
147 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
151template<
class TypeTag>
153{
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
156template<
class TypeTag>
174template <
class TypeTag>
197 wettingPhaseIdx = MaterialLaw::wettingPhaseIdx,
198 nonWettingPhaseIdx = MaterialLaw::nonWettingPhaseIdx,
201 numPhases = FluidSystem::numPhases,
204 dim = GridView::dimension,
205 dimWorld = GridView::dimensionworld
208 using FluidState = Opm::ImmiscibleFluidState<Scalar, FluidSystem>;
210 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
211 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
216 bool contains(Dune::GeometryType gt)
217 {
return gt.dim() == dim - 1; }
219 using FaceMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
228 : ParentType(simulator)
236 ParentType::finishInit();
239 temperature_ = 273.15 + 20;
241 matrixMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.0);
242 matrixMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
243 fractureMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.0);
244 fractureMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
247 matrixMaterialParams_.setEntryPC(0.0);
248 matrixMaterialParams_.setMaxPC(2000.0);
249 fractureMaterialParams_.setEntryPC(0.0);
250 fractureMaterialParams_.setMaxPC(1000.0);
254 matrixMaterialParams_.setEntryPressure(2000);
255 matrixMaterialParams_.setLambda(2.0);
256 matrixMaterialParams_.setPcLowSw(1e-1);
257 fractureMaterialParams_.setEntryPressure(1000);
258 fractureMaterialParams_.setLambda(2.0);
259 fractureMaterialParams_.setPcLowSw(5e-2);
263 matrixMaterialParams_.setVgAlpha(0.0037);
264 matrixMaterialParams_.setVgN(4.7);
265 fractureMaterialParams_.setVgAlpha(0.0025);
266 fractureMaterialParams_.setVgN(4.7);
269 matrixMaterialParams_.finalize();
270 fractureMaterialParams_.finalize();
272 matrixK_ = this->toDimMatrix_(1e-15);
273 fractureK_ = this->toDimMatrix_(1e5 * 1e-15);
275 matrixPorosity_ = 0.10;
276 fracturePorosity_ = 0.25;
277 fractureWidth_ = 1e-3;
280 initEnergyParams_(thermalConductionParams_, matrixPorosity_);
288 ParentType::registerParameters();
290 Parameters::SetDefault<Parameters::GridFile>(
"data/fracture.art.dgf");
291 Parameters::SetDefault<Parameters::EndTime<Scalar>>(3e3);
292 Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(100);
305 std::ostringstream oss;
306 oss <<
"fracture_" << Model::name();
322 this->model().globalStorage(storage);
325 if (this->gridView().comm().rank() == 0) {
326 std::cout <<
"Storage: " << storage << std::endl << std::flush;
334 template <
class Context>
336 [[maybe_unused]]
unsigned spaceIdx,
337 [[maybe_unused]]
unsigned timeIdx)
const
338 {
return temperature_; }
350 template <
class Context>
352 [[maybe_unused]]
unsigned spaceIdx,
353 [[maybe_unused]]
unsigned timeIdx)
const
361 template <
class Context>
363 [[maybe_unused]]
unsigned spaceIdx,
364 [[maybe_unused]]
unsigned timeIdx)
const
365 {
return fractureK_; }
370 template <
class Context>
371 Scalar
porosity([[maybe_unused]]
const Context& context,
372 [[maybe_unused]]
unsigned spaceIdx,
373 [[maybe_unused]]
unsigned timeIdx)
const
374 {
return matrixPorosity_; }
381 template <
class Context>
383 [[maybe_unused]]
unsigned spaceIdx,
384 [[maybe_unused]]
unsigned timeIdx)
const
385 {
return fracturePorosity_; }
390 template <
class Context>
392 [[maybe_unused]]
unsigned spaceIdx,
393 [[maybe_unused]]
unsigned timeIdx)
const
394 {
return matrixMaterialParams_; }
401 template <
class Context>
403 [[maybe_unused]]
unsigned spaceIdx,
404 [[maybe_unused]]
unsigned timeIdx)
const
405 {
return fractureMaterialParams_; }
411 {
return this->simulator().vanguard().fractureMapper(); }
425 template <
class Context>
427 [[maybe_unused]]
unsigned spaceIdx1,
428 [[maybe_unused]]
unsigned spaceIdx2,
429 [[maybe_unused]]
unsigned timeIdx)
const
430 {
return fractureWidth_; }
435 template <
class Context>
436 const ThermalConductionLawParams&
438 [[maybe_unused]]
unsigned spaceIdx,
439 [[maybe_unused]]
unsigned timeIdx)
const
440 {
return thermalConductionParams_; }
447 template <
class Context>
448 const SolidEnergyLawParams&
450 [[maybe_unused]]
unsigned spaceIdx,
451 [[maybe_unused]]
unsigned timeIdx)
const
452 {
return solidEnergyParams_; }
464 template <
class Context>
465 void boundary(BoundaryRateVector& values,
const Context& context,
466 unsigned spaceIdx,
unsigned timeIdx)
const
468 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
470 if (onRightBoundary_(pos)) {
473 FluidState fluidState;
474 fluidState.setTemperature(temperature_);
476 fluidState.setSaturation(wettingPhaseIdx, 0.0);
477 fluidState.setSaturation(nonWettingPhaseIdx,
478 1.0 - fluidState.saturation(wettingPhaseIdx));
480 fluidState.setPressure(wettingPhaseIdx, 1e5);
481 fluidState.setPressure(nonWettingPhaseIdx, fluidState.pressure(wettingPhaseIdx));
483 typename FluidSystem::template ParameterCache<Scalar> paramCache;
484 paramCache.updateAll(fluidState);
485 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
486 fluidState.setDensity(phaseIdx,
487 FluidSystem::density(fluidState, paramCache, phaseIdx));
488 fluidState.setViscosity(phaseIdx,
489 FluidSystem::viscosity(fluidState, paramCache, phaseIdx));
493 values.setFreeFlow(context, spaceIdx, timeIdx, fluidState);
511 template <
class Context>
513 unsigned spaceIdx,
unsigned timeIdx)
const
515 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
517 if (!onLeftBoundary_(pos))
521 unsigned globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
531 FluidState fractureFluidState;
532 fractureFluidState.setTemperature(temperature_ + 10.0);
534 fractureFluidState.setSaturation(wettingPhaseIdx, 1.0);
535 fractureFluidState.setSaturation(nonWettingPhaseIdx,
536 1.0 - fractureFluidState.saturation(
539 Scalar pCFracture[numPhases];
540 MaterialLaw::capillaryPressures(pCFracture, fractureMaterialParams_,
543 fractureFluidState.setPressure(wettingPhaseIdx, 1.0e5);
544 fractureFluidState.setPressure(nonWettingPhaseIdx,
545 fractureFluidState.pressure(wettingPhaseIdx)
546 + (pCFracture[nonWettingPhaseIdx]
547 - pCFracture[wettingPhaseIdx]));
550 constraints.assignNaiveFromFracture(fractureFluidState,
551 matrixMaterialParams_);
557 template <
class Context>
559 [[maybe_unused]]
const Context& context,
560 [[maybe_unused]]
unsigned spaceIdx,
561 [[maybe_unused]]
unsigned timeIdx)
const
563 FluidState fluidState;
564 fluidState.setTemperature(temperature_);
565 fluidState.setPressure(FluidSystem::wettingPhaseIdx, 1e5);
566 fluidState.setPressure(nonWettingPhaseIdx, fluidState.pressure(wettingPhaseIdx));
568 fluidState.setSaturation(wettingPhaseIdx, 0.0);
569 fluidState.setSaturation(nonWettingPhaseIdx,
570 1.0 - fluidState.saturation(wettingPhaseIdx));
572 values.assignNaive(fluidState);
581 template <
class Context>
583 [[maybe_unused]]
const Context& context,
584 [[maybe_unused]]
unsigned spaceIdx,
585 [[maybe_unused]]
unsigned timeIdx)
const
586 { rate = Scalar(0.0); }
591 bool onLeftBoundary_(
const GlobalPosition& pos)
const
592 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
594 bool onRightBoundary_(
const GlobalPosition& pos)
const
595 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
597 bool onLowerBoundary_(
const GlobalPosition& pos)
const
598 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
600 bool onUpperBoundary_(
const GlobalPosition& pos)
const
601 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
603 void initEnergyParams_(ThermalConductionLawParams& params, Scalar poro)
606 solidEnergyParams_.setSolidHeatCapacity(790.0
608 solidEnergyParams_.finalize();
610 Scalar lambdaGranite = 2.8;
613 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
614 fs.setTemperature(293.15);
615 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
616 fs.setPressure(phaseIdx, 1.0135e5);
619 typename FluidSystem::template ParameterCache<Scalar> paramCache;
620 paramCache.updateAll(fs);
621 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
622 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
623 fs.setDensity(phaseIdx, rho);
626 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
627 Scalar lambdaSaturated;
628 if (FluidSystem::isLiquid(phaseIdx)) {
629 Scalar lambdaFluid = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
631 std::pow(lambdaGranite, (1 - poro))
632 + std::pow(lambdaFluid, poro);
635 lambdaSaturated = std::pow(lambdaGranite, (1 - poro));
637 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
640 Scalar lambdaVac = std::pow(lambdaGranite, (1 - poro));
641 params.setVacuumLambda(lambdaVac);
645 DimMatrix fractureK_;
647 Scalar matrixPorosity_;
648 Scalar fracturePorosity_;
650 Scalar fractureWidth_;
652 MaterialLawParams fractureMaterialParams_;
653 MaterialLawParams matrixMaterialParams_;
655 ThermalConductionLawParams thermalConductionParams_;
656 SolidEnergyLawParams solidEnergyParams_;
Provides a simulator vanguard which creates a grid by parsing a Dune Grid Format (DGF) file.
Definition: dgfvanguard.hh:50
Stores the topology of fractures.
Definition: fracturemapper.hh:43
Two-phase problem which involves fractures.
Definition: fractureproblem.hh:176
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:449
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:351
Scalar fractureWidth(const Context &context, unsigned spaceIdx1, unsigned spaceIdx2, unsigned timeIdx) const
Returns the width of the fracture.
Definition: fractureproblem.hh:426
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: fractureproblem.hh:558
static void registerParameters()
Definition: fractureproblem.hh:286
Scalar fracturePorosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
The porosity inside the fractures.
Definition: fractureproblem.hh:382
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: fractureproblem.hh:234
std::string name() const
The problem name.
Definition: fractureproblem.hh:303
FractureProblem(Simulator &simulator)
Definition: fractureproblem.hh:227
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the constraints for a control volume.
Definition: fractureproblem.hh:512
const DimMatrix & fractureIntrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Intrinsic permeability of fractures.
Definition: fractureproblem.hh:362
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:391
void endTimeStep()
Called directly after the time integration.
Definition: fractureproblem.hh:313
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:437
void source(RateVector &rate, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: fractureproblem.hh:582
const MaterialLawParams & fractureMaterialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
The parameters for the material law inside the fractures.
Definition: fractureproblem.hh:402
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:335
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: fractureproblem.hh:465
const FractureMapper & fractureMapper() const
Returns the object representating the fracture topology.
Definition: fractureproblem.hh:410
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fractureproblem.hh:371
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
@ simplex
Definition: vcfvstencil.hh:54
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
Specify whether the some degrees of fredom can be constraint.
Definition: fvbaseproperties.hh:199
Specify whether energy should be considered as a conservation quantity or not.
Definition: multiphasebaseproperties.hh:76
Dune::ALUGrid< 2, 2, Dune::simplex, Dune::nonconforming > type
Definition: fractureproblem.hh:78
The type of the DUNE grid.
Definition: basicproperties.hh:100
Opm::EffToAbsLaw< EffectiveLaw > type
Definition: fractureproblem.hh:130
The material law which ought to be used (extracted from the spatial parameters)
Definition: multiphasebaseproperties.hh:51
Opm::LiquidPhase< Scalar, Opm::DNAPL< Scalar > > type
Definition: fractureproblem.hh:107
The non-wetting phase for two-phase models.
Definition: immiscibleproperties.hh:44
The type of the problem.
Definition: fvbaseproperties.hh:81
Opm::ConstantSolidHeatCapLaw< GetPropType< TypeTag, Properties::Scalar > > type
Definition: fractureproblem.hh:153
The material law for the energy stored in the solid matrix.
Definition: multiphasebaseproperties.hh:57
Definition: fractureproblem.hh:72
std::tuple< DiscreteFractureModel > InheritsFrom
Definition: fractureproblem.hh:72
Opm::SomertonThermalConductionLaw< FluidSystem, Scalar > type
Definition: fractureproblem.hh:147
The material law for thermal conduction.
Definition: multiphasebaseproperties.hh:63
Property which provides a Vanguard (manages grids)
Definition: basicproperties.hh:96
Opm::LiquidPhase< Scalar, Opm::SimpleH2O< Scalar > > type
Definition: fractureproblem.hh:96
The wetting phase for two-phase models.
Definition: immiscibleproperties.hh:41