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/heatconduction/Somerton.hpp>
49 #include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp>
50 #include <opm/material/components/SimpleH2O.hpp>
51 #include <opm/material/components/Dnapl.hpp>
52 #include <opm/common/Unused.hpp>
54 #include <dune/common/version.hh>
55 #include <dune/common/fmatrix.hh>
56 #include <dune/common/fvector.hh>
63 template <
class TypeTag>
68 namespace Properties {
75 Dune::ALUGrid</*dim=*/2, /*dimWorld=*/2, Dune::simplex, Dune::nonconforming>);
90 typedef Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> > type;
100 typedef Opm::LiquidPhase<Scalar, Opm::DNAPL<Scalar> > type;
104 SET_PROP(FractureProblem, MaterialLaw)
107 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
108 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
109 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
112 typedef Opm::TwoPhaseMaterialTraits<Scalar,
113 FluidSystem::wettingPhaseIdx,
114 FluidSystem::nonWettingPhaseIdx>
119 typedef Opm::RegularizedBrooksCorey<Traits> EffectiveLaw;
123 typedef Opm::EffToAbsLaw<EffectiveLaw> type;
130 SET_PROP(FractureProblem, HeatConductionLaw)
134 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
138 typedef Opm::Somerton<FluidSystem, Scalar> type;
171 template <
class TypeTag>
172 class FractureProblem :
public GET_PROP_TYPE(TypeTag, BaseProblem)
174 typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
176 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
177 typedef typename
GET_PROP_TYPE(TypeTag, WettingPhase) WettingPhase;
178 typedef typename
GET_PROP_TYPE(TypeTag, NonwettingPhase) NonwettingPhase;
179 typedef typename
GET_PROP_TYPE(TypeTag, Constraints) Constraints;
181 typedef typename
GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
182 typedef typename
GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
183 typedef typename
GET_PROP_TYPE(TypeTag, RateVector) RateVector;
184 typedef typename
GET_PROP_TYPE(TypeTag, Simulator) Simulator;
186 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
187 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
188 typedef typename
GET_PROP_TYPE(TypeTag, HeatConductionLawParams) HeatConductionLawParams;
193 wettingPhaseIdx = MaterialLaw::wettingPhaseIdx,
194 nonWettingPhaseIdx = MaterialLaw::nonWettingPhaseIdx,
197 numPhases = FluidSystem::numPhases,
200 dim = GridView::dimension,
201 dimWorld = GridView::dimensionworld
204 typedef Opm::ImmiscibleFluidState<Scalar, FluidSystem> FluidState;
206 typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
207 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
212 bool contains(Dune::GeometryType gt)
213 {
return gt.dim() == dim - 1; }
215 typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, FaceLayout> FaceMapper;
224 : ParentType(simulator)
232 ParentType::finishInit();
235 temperature_ = 273.15 + 20;
237 matrixMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.0);
238 matrixMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
239 fractureMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.0);
240 fractureMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
243 matrixMaterialParams_.setEntryPC(0.0);
244 matrixMaterialParams_.setMaxPC(2000.0);
245 fractureMaterialParams_.setEntryPC(0.0);
246 fractureMaterialParams_.setMaxPC(1000.0);
249 #if 1 // Brooks-Corey
250 matrixMaterialParams_.setEntryPressure(2000);
251 matrixMaterialParams_.setLambda(2.0);
252 matrixMaterialParams_.setPcLowSw(1e-1);
253 fractureMaterialParams_.setEntryPressure(1000);
254 fractureMaterialParams_.setLambda(2.0);
255 fractureMaterialParams_.setPcLowSw(5e-2);
258 #if 0 // van Genuchten
259 matrixMaterialParams_.setVgAlpha(0.0037);
260 matrixMaterialParams_.setVgN(4.7);
261 fractureMaterialParams_.setVgAlpha(0.0025);
262 fractureMaterialParams_.setVgN(4.7);
265 matrixMaterialParams_.finalize();
266 fractureMaterialParams_.finalize();
268 matrixK_ = this->toDimMatrix_(1e-15);
269 fractureK_ = this->toDimMatrix_(1e5 * 1e-15);
271 matrixPorosity_ = 0.10;
272 fracturePorosity_ = 0.25;
273 fractureWidth_ = 1e-3;
276 computeHeatCondParams_(heatCondParams_, matrixPorosity_);
289 std::ostringstream oss;
290 oss <<
"fracture_" << Model::name();
306 this->model().globalStorage(storage);
309 if (this->gridView().comm().rank() == 0) {
310 std::cout <<
"Storage: " << storage << std::endl << std::flush;
318 template <
class Context>
320 unsigned spaceIdx OPM_UNUSED,
321 unsigned timeIdx OPM_UNUSED)
const
322 {
return temperature_; }
334 template <
class Context>
336 unsigned spaceIdx OPM_UNUSED,
337 unsigned timeIdx OPM_UNUSED)
const
345 template <
class Context>
347 unsigned spaceIdx OPM_UNUSED,
348 unsigned timeIdx OPM_UNUSED)
const
349 {
return fractureK_; }
354 template <
class Context>
356 unsigned spaceIdx OPM_UNUSED,
357 unsigned timeIdx OPM_UNUSED)
const
358 {
return matrixPorosity_; }
365 template <
class Context>
367 unsigned spaceIdx OPM_UNUSED,
368 unsigned timeIdx OPM_UNUSED)
const
369 {
return fracturePorosity_; }
374 template <
class Context>
376 unsigned spaceIdx OPM_UNUSED,
377 unsigned timeIdx OPM_UNUSED)
const
378 {
return matrixMaterialParams_; }
385 template <
class Context>
387 unsigned spaceIdx OPM_UNUSED,
388 unsigned timeIdx OPM_UNUSED)
const
389 {
return fractureMaterialParams_; }
395 {
return this->simulator().gridManager().fractureMapper(); }
409 template <
class Context>
411 unsigned spaceIdx1 OPM_UNUSED,
412 unsigned spaceIdx2 OPM_UNUSED,
413 unsigned timeIdx OPM_UNUSED)
const
414 {
return fractureWidth_; }
419 template <
class Context>
420 const HeatConductionLawParams &
422 unsigned spaceIdx OPM_UNUSED,
423 unsigned timeIdx OPM_UNUSED)
const
424 {
return heatCondParams_; }
431 template <
class Context>
433 unsigned spaceIdx OPM_UNUSED,
434 unsigned timeIdx OPM_UNUSED)
const
450 template <
class Context>
451 void boundary(BoundaryRateVector& values,
const Context& context,
452 unsigned spaceIdx,
unsigned timeIdx)
const
454 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
456 if (onRightBoundary_(pos)) {
459 FluidState fluidState;
460 fluidState.setTemperature(temperature_);
462 fluidState.setSaturation(wettingPhaseIdx, 0.0);
463 fluidState.setSaturation(nonWettingPhaseIdx,
464 1.0 - fluidState.saturation(wettingPhaseIdx));
466 fluidState.setPressure(wettingPhaseIdx, 1e5);
467 fluidState.setPressure(nonWettingPhaseIdx, fluidState.pressure(wettingPhaseIdx));
470 values.setFreeFlow(context, spaceIdx, timeIdx, fluidState);
488 template <
class Context>
489 void constraints(Constraints& constraints,
const Context& context,
490 unsigned spaceIdx,
unsigned timeIdx)
const
492 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
494 if (!onLeftBoundary_(pos))
498 unsigned globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
499 if (!fractureMapper().isFractureVertex(globalIdx)) {
508 FluidState fractureFluidState;
509 fractureFluidState.setTemperature(temperature_ + 10.0);
511 fractureFluidState.setSaturation(wettingPhaseIdx, 1.0);
512 fractureFluidState.setSaturation(nonWettingPhaseIdx,
513 1.0 - fractureFluidState.saturation(
516 Scalar pCFracture[numPhases];
517 MaterialLaw::capillaryPressures(pCFracture, fractureMaterialParams_,
520 fractureFluidState.setPressure(wettingPhaseIdx, 1.0e5);
521 fractureFluidState.setPressure(nonWettingPhaseIdx,
522 fractureFluidState.pressure(wettingPhaseIdx)
523 + (pCFracture[nonWettingPhaseIdx]
524 - pCFracture[wettingPhaseIdx]));
526 constraints.setActive(
true);
527 constraints.assignNaiveFromFracture(fractureFluidState,
528 matrixMaterialParams_);
534 template <
class Context>
536 const Context& context OPM_UNUSED,
537 unsigned spaceIdx OPM_UNUSED,
538 unsigned timeIdx OPM_UNUSED)
const
540 FluidState fluidState;
541 fluidState.setTemperature(temperature_);
542 fluidState.setPressure(FluidSystem::wettingPhaseIdx, 1e5);
543 fluidState.setPressure(nonWettingPhaseIdx, fluidState.pressure(wettingPhaseIdx));
545 fluidState.setSaturation(wettingPhaseIdx, 0.0);
546 fluidState.setSaturation(nonWettingPhaseIdx,
547 1.0 - fluidState.saturation(wettingPhaseIdx));
549 values.assignNaive(fluidState);
558 template <
class Context>
560 const Context& context OPM_UNUSED,
561 unsigned spaceIdx OPM_UNUSED,
562 unsigned timeIdx OPM_UNUSED)
const
563 { rate = Scalar(0.0); }
568 bool onLeftBoundary_(
const GlobalPosition& pos)
const
569 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
571 bool onRightBoundary_(
const GlobalPosition& pos)
const
572 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
574 bool onLowerBoundary_(
const GlobalPosition& pos)
const
575 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
577 bool onUpperBoundary_(
const GlobalPosition& pos)
const
578 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
580 void computeHeatCondParams_(HeatConductionLawParams& params, Scalar poro)
582 Scalar lambdaGranite = 2.8;
585 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
586 fs.setTemperature(293.15);
587 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
588 fs.setPressure(phaseIdx, 1.0135e5);
591 typename FluidSystem::template ParameterCache<Scalar> paramCache;
592 paramCache.updateAll(fs);
593 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
594 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
595 fs.setDensity(phaseIdx, rho);
598 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
599 Scalar lambdaSaturated;
600 if (FluidSystem::isLiquid(phaseIdx)) {
601 Scalar lambdaFluid = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
603 std::pow(lambdaGranite, (1 - poro))
604 + std::pow(lambdaFluid, poro);
607 lambdaSaturated = std::pow(lambdaGranite, (1 - poro));
609 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
612 Scalar lambdaVac = std::pow(lambdaGranite, (1 - poro));
613 params.setVacuumLambda(lambdaVac);
617 DimMatrix fractureK_;
619 Scalar matrixPorosity_;
620 Scalar fracturePorosity_;
622 Scalar fractureWidth_;
624 MaterialLawParams fractureMaterialParams_;
625 MaterialLawParams matrixMaterialParams_;
627 HeatConductionLawParams heatCondParams_;
634 #endif // EWOMS_FRACTURE_PROBLEM_HH
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: fractureproblem.hh:451
Stores the topology of fractures.
Definition: fracturemapper.hh:42
Scalar fracturePorosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
The porosity inside the fractures.
Definition: fractureproblem.hh:366
Two-phase problem which involves fractures.
Definition: fractureproblem.hh:64
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the constraints for a control volume.
Definition: fractureproblem.hh:489
A fully-implicit multi-phase flow model which assumes immiscibility of the phases and is able to incl...
Definition: discretefracturemodel.hh:50
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: fractureproblem.hh:559
void endTimeStep()
Called directly after the time integration.
Definition: fractureproblem.hh:297
const HeatConductionLawParams & heatConductionParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:421
void initial(PrimaryVariables &values, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the initial value for a control volume.
Definition: fractureproblem.hh:535
#define INHERITS_FROM(...)
Syntactic sugar for NEW_TYPE_TAG.
Definition: propertysystem.hh:230
FractureProblem(Simulator &simulator)
Definition: fractureproblem.hh:223
Scalar heatCapacitySolid(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:432
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:319
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: fractureproblem.hh:230
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:355
std::string name() const
The problem name.
Definition: fractureproblem.hh:287
Provides a grid manager which reads Dune Grid Format (DGF) files.
Definition: dgfgridmanager.hh:56
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
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:386
const DimMatrix & fractureIntrinsicPermeability(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Intrinsic permeability of fractures.
Definition: fractureproblem.hh:346
const DimMatrix & intrinsicPermeability(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:335
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
#define SET_STRING_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant string value.
Definition: propertysystem.hh:416
A fully-implicit multi-phase flow model which assumes immiscibility of the phases and is able to incl...
const MaterialLawParams & materialLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:375
Provides a grid manager which reads Dune Grid Format (DGF) files.
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:410
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394
const FractureMapper & fractureMapper() const
Returns the object representating the fracture topology.
Definition: fractureproblem.hh:394