My Project
obstacleproblem.hh
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
28 #ifndef EWOMS_OBSTACLE_PROBLEM_HH
29 #define EWOMS_OBSTACLE_PROBLEM_HH
30 
31 #include <opm/models/ncp/ncpproperties.hh>
32 
33 #include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
34 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
35 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
36 #include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
37 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
38 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
39 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
40 #include <opm/material/thermal/ConstantSolidHeatCapLaw.hpp>
41 #include <opm/material/thermal/SomertonThermalConductionLaw.hpp>
42 
43 #include <dune/grid/yaspgrid.hh>
44 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
45 
46 #include <dune/common/version.hh>
47 #include <dune/common/fvector.hh>
48 #include <dune/common/fmatrix.hh>
49 
50 #include <sstream>
51 #include <string>
52 #include <iostream>
53 
54 namespace Opm {
55 template <class TypeTag>
56 class ObstacleProblem;
57 }
58 
59 namespace Opm::Properties {
60 
61 namespace TTag {
63 }
64 
65 // Set the grid type
66 template<class TypeTag>
67 struct Grid<TypeTag, TTag::ObstacleBaseProblem> { using type = Dune::YaspGrid<2>; };
68 
69 // Set the problem property
70 template<class TypeTag>
71 struct Problem<TypeTag, TTag::ObstacleBaseProblem> { using type = Opm::ObstacleProblem<TypeTag>; };
72 
73 // Set fluid configuration
74 template<class TypeTag>
75 struct FluidSystem<TypeTag, TTag::ObstacleBaseProblem>
76 { using type = Opm::H2ON2FluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
77 
78 // Set the material Law
79 template<class TypeTag>
80 struct MaterialLaw<TypeTag, TTag::ObstacleBaseProblem>
81 {
82 private:
83  // define the material law
84  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
85  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
86  using MaterialTraits = Opm::TwoPhaseMaterialTraits<Scalar,
87  /*wettingPhaseIdx=*/FluidSystem::liquidPhaseIdx,
88  /*nonWettingPhaseIdx=*/FluidSystem::gasPhaseIdx>;
89 
90  using EffMaterialLaw = Opm::LinearMaterial<MaterialTraits>;
91 
92 public:
93  using type = Opm::EffToAbsLaw<EffMaterialLaw>;
94 };
95 
96 // Set the thermal conduction law
97 template<class TypeTag>
98 struct ThermalConductionLaw<TypeTag, TTag::ObstacleBaseProblem>
99 {
100 private:
101  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
102  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
103 
104 public:
105  // define the material law parameterized by absolute saturations
106  using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
107 };
108 
109 // set the energy storage law for the solid phase
110 template<class TypeTag>
111 struct SolidEnergyLaw<TypeTag, TTag::ObstacleBaseProblem>
112 { using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
113 
114 // Enable gravity
115 template<class TypeTag>
116 struct EnableGravity<TypeTag, TTag::ObstacleBaseProblem> { static constexpr bool value = true; };
117 
118 // The default for the end time of the simulation
119 template<class TypeTag>
120 struct EndTime<TypeTag, TTag::ObstacleBaseProblem>
121 {
122  using type = GetPropType<TypeTag, Scalar>;
123  static constexpr type value = 1e4;
124 };
125 
126 // The default for the initial time step size of the simulation
127 template<class TypeTag>
128 struct InitialTimeStepSize<TypeTag, TTag::ObstacleBaseProblem>
129 {
130  using type = GetPropType<TypeTag, Scalar>;
131  static constexpr type value = 250;
132 };
133 
134 // The default DGF file to load
135 template<class TypeTag>
136 struct GridFile<TypeTag, TTag::ObstacleBaseProblem> { static constexpr auto value = "./data/obstacle_24x16.dgf"; };
137 
138 } // namespace Opm::Properties
139 
140 namespace Opm {
167 template <class TypeTag>
168 class ObstacleProblem : public GetPropType<TypeTag, Properties::BaseProblem>
169 {
170  using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
171 
172  using GridView = GetPropType<TypeTag, Properties::GridView>;
173  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
174  using EqVector = GetPropType<TypeTag, Properties::EqVector>;
175  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
176  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
177  using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
178  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
179  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
180  using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
181  using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
182  using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
183 
184  enum {
185  // Grid and world dimension
186  dim = GridView::dimension,
187  dimWorld = GridView::dimensionworld,
188  numPhases = getPropValue<TypeTag, Properties::NumPhases>(),
189  gasPhaseIdx = FluidSystem::gasPhaseIdx,
190  liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
191  H2OIdx = FluidSystem::H2OIdx,
192  N2Idx = FluidSystem::N2Idx
193  };
194 
195  using GlobalPosition = Dune::FieldVector<typename GridView::ctype, dimWorld>;
196  using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
197  using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
198  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
199  using Model = GetPropType<TypeTag, Properties::Model>;
200 
201 public:
205  ObstacleProblem(Simulator& simulator)
206  : ParentType(simulator)
207  { }
208 
212  void finishInit()
213  {
214  ParentType::finishInit();
215 
216  eps_ = 1e-6;
217  temperature_ = 273.15 + 25; // -> 25°C
218 
219  // initialize the tables of the fluid system
220  Scalar Tmin = temperature_ - 1.0;
221  Scalar Tmax = temperature_ + 1.0;
222  unsigned nT = 3;
223 
224  Scalar pmin = 1.0e5 * 0.75;
225  Scalar pmax = 2.0e5 * 1.25;
226  unsigned np = 1000;
227 
228  FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
229 
230  // intrinsic permeabilities
231  coarseK_ = this->toDimMatrix_(1e-12);
232  fineK_ = this->toDimMatrix_(1e-15);
233 
234  // the porosity
235  finePorosity_ = 0.3;
236  coarsePorosity_ = 0.3;
237 
238  // residual saturations
239  fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
240  fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
241  coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
242  coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
243 
244  // parameters for the linear law, i.e. minimum and maximum
245  // pressures
246  fineMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
247  fineMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
248  coarseMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
249  coarseMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
250 
251  /*
252  // entry pressures for Brooks-Corey
253  fineMaterialParams_.setEntryPressure(5e3);
254  coarseMaterialParams_.setEntryPressure(1e3);
255 
256  // Brooks-Corey shape parameters
257  fineMaterialParams_.setLambda(2);
258  coarseMaterialParams_.setLambda(2);
259  */
260 
261  fineMaterialParams_.finalize();
262  coarseMaterialParams_.finalize();
263 
264  // parameters for the somerton law of thermal conduction
265  computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
266  computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
267 
268  // assume constant volumetric heat capacity and granite
269  solidEnergyLawParams_.setSolidHeatCapacity(790.0 // specific heat capacity of granite [J / (kg K)]
270  * 2700.0); // density of granite [kg/m^3]
271  solidEnergyLawParams_.finalize();
272 
273  initFluidStates_();
274  }
275 
279  void endTimeStep()
280  {
281 #ifndef NDEBUG
282  this->model().checkConservativeness();
283 
284  // Calculate storage terms of the individual phases
285  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
286  PrimaryVariables phaseStorage;
287  this->model().globalPhaseStorage(phaseStorage, phaseIdx);
288 
289  if (this->gridView().comm().rank() == 0) {
290  std::cout << "Storage in " << FluidSystem::phaseName(phaseIdx)
291  << "Phase: [" << phaseStorage << "]"
292  << "\n" << std::flush;
293  }
294  }
295 
296  // Calculate total storage terms
297  EqVector storage;
298  this->model().globalStorage(storage);
299 
300  // Write mass balance information for rank 0
301  if (this->gridView().comm().rank() == 0) {
302  std::cout << "Storage total: [" << storage << "]"
303  << "\n" << std::flush;
304  }
305 #endif // NDEBUG
306  }
307 
312 
316  std::string name() const
317  {
318  std::ostringstream oss;
319  oss << "obstacle"
320  << "_" << Model::name();
321  return oss.str();
322  }
323 
329  template <class Context>
330  Scalar temperature(const Context& /*context*/,
331  unsigned /*spaceIdx*/,
332  unsigned /*timeIdx*/) const
333  { return temperature_; }
334 
338  template <class Context>
339  const DimMatrix&
340  intrinsicPermeability(const Context& context,
341  unsigned spaceIdx,
342  unsigned timeIdx) const
343  {
344  if (isFineMaterial_(context.pos(spaceIdx, timeIdx)))
345  return fineK_;
346  return coarseK_;
347  }
348 
352  template <class Context>
353  Scalar porosity(const Context& context,
354  unsigned spaceIdx,
355  unsigned timeIdx) const
356  {
357  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
358  if (isFineMaterial_(pos))
359  return finePorosity_;
360  else
361  return coarsePorosity_;
362  }
363 
367  template <class Context>
368  const MaterialLawParams&
369  materialLawParams(const Context& context,
370  unsigned spaceIdx,
371  unsigned timeIdx) const
372  {
373  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
374  if (isFineMaterial_(pos))
375  return fineMaterialParams_;
376  else
377  return coarseMaterialParams_;
378  }
379 
385  template <class Context>
386  const SolidEnergyLawParams&
387  solidEnergyLawParams(const Context& /*context*/,
388  unsigned /*spaceIdx*/,
389  unsigned /*timeIdx*/) const
390  { return solidEnergyLawParams_; }
391 
395  template <class Context>
396  const ThermalConductionLawParams &
397  thermalConductionParams(const Context& context,
398  unsigned spaceIdx,
399  unsigned timeIdx) const
400  {
401  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
402  if (isFineMaterial_(pos))
403  return fineThermalCondParams_;
404  return coarseThermalCondParams_;
405  }
406 
408 
413 
417  template <class Context>
418  void boundary(BoundaryRateVector& values,
419  const Context& context,
420  unsigned spaceIdx,
421  unsigned timeIdx) const
422  {
423  const auto& pos = context.pos(spaceIdx, timeIdx);
424 
425  if (onInlet_(pos))
426  values.setFreeFlow(context, spaceIdx, timeIdx, inletFluidState_);
427  else if (onOutlet_(pos))
428  values.setFreeFlow(context, spaceIdx, timeIdx, outletFluidState_);
429  else
430  values.setNoFlow();
431  }
432 
434 
439 
443  template <class Context>
444  void initial(PrimaryVariables& values,
445  const Context& context,
446  unsigned spaceIdx,
447  unsigned timeIdx) const
448  {
449  const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
450  values.assignMassConservative(outletFluidState_, matParams);
451  }
452 
459  template <class Context>
460  void source(RateVector& rate,
461  const Context& /*context*/,
462  unsigned /*spaceIdx*/,
463  unsigned /*timeIdx*/) const
464  { rate = 0.0; }
465 
467 
468 private:
473  bool isFineMaterial_(const GlobalPosition& pos) const
474  { return 10 <= pos[0] && pos[0] <= 20 && 0 <= pos[1] && pos[1] <= 35; }
475 
476  bool onInlet_(const GlobalPosition& globalPos) const
477  {
478  Scalar x = globalPos[0];
479  Scalar y = globalPos[1];
480  return x >= 60 - eps_ && y <= 10;
481  }
482 
483  bool onOutlet_(const GlobalPosition& globalPos) const
484  {
485  Scalar x = globalPos[0];
486  Scalar y = globalPos[1];
487  return x < eps_ && y <= 10;
488  }
489 
490  void initFluidStates_()
491  {
492  initFluidState_(inletFluidState_, coarseMaterialParams_,
493  /*isInlet=*/true);
494  initFluidState_(outletFluidState_, coarseMaterialParams_,
495  /*isInlet=*/false);
496  }
497 
498  template <class FluidState>
499  void initFluidState_(FluidState& fs, const MaterialLawParams& matParams, bool isInlet)
500  {
501  unsigned refPhaseIdx;
502  unsigned otherPhaseIdx;
503 
504  // set the fluid temperatures
505  fs.setTemperature(temperature_);
506 
507  if (isInlet) {
508  // only liquid on inlet
509  refPhaseIdx = liquidPhaseIdx;
510  otherPhaseIdx = gasPhaseIdx;
511 
512  // set liquid saturation
513  fs.setSaturation(liquidPhaseIdx, 1.0);
514 
515  // set pressure of the liquid phase
516  fs.setPressure(liquidPhaseIdx, 2e5);
517 
518  // set the liquid composition to pure water
519  fs.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
520  fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
521  }
522  else {
523  // elsewhere, only gas
524  refPhaseIdx = gasPhaseIdx;
525  otherPhaseIdx = liquidPhaseIdx;
526 
527  // set gas saturation
528  fs.setSaturation(gasPhaseIdx, 1.0);
529 
530  // set pressure of the gas phase
531  fs.setPressure(gasPhaseIdx, 1e5);
532 
533  // set the gas composition to 99% nitrogen and 1% steam
534  fs.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
535  fs.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
536  }
537 
538  // set the other saturation
539  fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
540 
541  // calulate the capillary pressure
542  PhaseVector pC;
543  MaterialLaw::capillaryPressures(pC, matParams, fs);
544  fs.setPressure(otherPhaseIdx, fs.pressure(refPhaseIdx)
545  + (pC[otherPhaseIdx] - pC[refPhaseIdx]));
546 
547  // make the fluid state consistent with local thermodynamic
548  // equilibrium
549  using ComputeFromReferencePhase = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
550 
551  typename FluidSystem::template ParameterCache<Scalar> paramCache;
552  ComputeFromReferencePhase::solve(fs, paramCache, refPhaseIdx,
553  /*setViscosity=*/true,
554  /*setEnthalpy=*/false);
555  }
556 
557  void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
558  {
559  Scalar lambdaWater = 0.6;
560  Scalar lambdaGranite = 2.8;
561 
562  Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
563  * std::pow(lambdaWater, poro);
564  Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
565 
566  params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
567  params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
568  params.setVacuumLambda(lambdaDry);
569  }
570 
571  DimMatrix coarseK_;
572  DimMatrix fineK_;
573 
574  Scalar coarsePorosity_;
575  Scalar finePorosity_;
576 
577  MaterialLawParams fineMaterialParams_;
578  MaterialLawParams coarseMaterialParams_;
579 
580  ThermalConductionLawParams fineThermalCondParams_;
581  ThermalConductionLawParams coarseThermalCondParams_;
582  SolidEnergyLawParams solidEnergyLawParams_;
583 
584  Opm::CompositionalFluidState<Scalar, FluidSystem> inletFluidState_;
585  Opm::CompositionalFluidState<Scalar, FluidSystem> outletFluidState_;
586 
587  Scalar temperature_;
588  Scalar eps_;
589 };
590 } // namespace Opm
591 
592 #endif
Problem where liquid water is first stopped by a low-permeability lens and then seeps though it.
Definition: obstacleproblem.hh:169
const SolidEnergyLawParams & solidEnergyLawParams(const Context &, unsigned, unsigned) const
Return the parameters for the energy storage law of the rock.
Definition: obstacleproblem.hh:387
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:353
void finishInit()
Definition: obstacleproblem.hh:212
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:369
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: obstacleproblem.hh:330
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:340
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:418
std::string name() const
Definition: obstacleproblem.hh:316
void endTimeStep()
Definition: obstacleproblem.hh:279
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:444
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: obstacleproblem.hh:460
ObstacleProblem(Simulator &simulator)
Definition: obstacleproblem.hh:205
const ThermalConductionLawParams & thermalConductionParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:397
Definition: obstacleproblem.hh:62