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;
94 SET_PROP(FractureProblem, NonwettingPhase)
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;
175 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
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;
180 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
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;
185 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
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;
189 typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
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>
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);
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]));
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 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
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
const HeatConductionLawParams & heatConductionParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:421
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: fractureproblem.hh:451
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
Definition: baseauxiliarymodule.hh:37
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:355
Stores the topology of fractures.
Definition: fracturemapper.hh:42
Two-phase problem which involves fractures.
Definition: fractureproblem.hh:64
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
Scalar heatCapacitySolid(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:432
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
void endTimeStep()
Called directly after the time integration.
Definition: fractureproblem.hh:297
#define INHERITS_FROM(...)
Syntactic sugar for NEW_TYPE_TAG.
Definition: propertysystem.hh:230
FractureProblem(Simulator &simulator)
Definition: fractureproblem.hh:223
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the constraints for a control volume.
Definition: fractureproblem.hh:489
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: fractureproblem.hh:230
Provides a grid manager which reads Dune Grid Format (DGF) files.
Definition: dgfgridmanager.hh:56
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
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:319
const FractureMapper & fractureMapper() const
Returns the object representating the fracture topology.
Definition: fractureproblem.hh:394
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
std::string name() const
The problem name.
Definition: fractureproblem.hh:287
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 MaterialLawParams & materialLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:375
const DimMatrix & intrinsicPermeability(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: fractureproblem.hh:335
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
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
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...
Provides a grid manager which reads Dune Grid Format (DGF) files.
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394