waterairproblem.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_WATER_AIR_PROBLEM_HH
29 #define EWOMS_WATER_AIR_PROBLEM_HH
30 
33 
34 #include <opm/material/fluidsystems/H2OAirFluidSystem.hpp>
35 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
36 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
37 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
38 #include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
39 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
40 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
41 #include <opm/material/heatconduction/Somerton.hpp>
42 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
43 #include <opm/common/Unused.hpp>
44 
45 #include <dune/grid/yaspgrid.hh>
46 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
47 
48 #include <dune/common/fvector.hh>
49 #include <dune/common/fmatrix.hh>
50 #include <dune/common/version.hh>
51 
52 #include <sstream>
53 #include <string>
54 
55 namespace Ewoms {
56 template <class TypeTag>
58 }
59 
60 namespace Ewoms {
61 namespace Properties {
62 NEW_TYPE_TAG(WaterAirBaseProblem);
63 
64 // Set the grid type
65 SET_TYPE_PROP(WaterAirBaseProblem, Grid, Dune::YaspGrid<2>);
66 
67 // Set the problem property
68 SET_TYPE_PROP(WaterAirBaseProblem, Problem, Ewoms::WaterAirProblem<TypeTag>);
69 
70 // Set the material Law
71 SET_PROP(WaterAirBaseProblem, MaterialLaw)
72 {
73 private:
74  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
75  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
76  typedef Opm::TwoPhaseMaterialTraits<Scalar,
77  /*wettingPhaseIdx=*/FluidSystem::liquidPhaseIdx,
78  /*nonWettingPhaseIdx=*/FluidSystem::gasPhaseIdx> Traits;
79 
80  // define the material law which is parameterized by effective
81  // saturations
82  typedef Opm::RegularizedBrooksCorey<Traits> EffMaterialLaw;
83 
84 public:
85  // define the material law parameterized by absolute saturations
86  // which uses the two-phase API
87  typedef Opm::EffToAbsLaw<EffMaterialLaw> type;
88 };
89 
90 // Set the heat conduction law
91 SET_PROP(WaterAirBaseProblem, HeatConductionLaw)
92 {
93 private:
94  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
95  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
96 
97 public:
98  // define the material law parameterized by absolute saturations
99  typedef Opm::Somerton<FluidSystem, Scalar> type;
100 };
101 
102 // Set the fluid system. in this case, we use the one which describes
103 // air and water
104 SET_TYPE_PROP(WaterAirBaseProblem, FluidSystem,
105  Opm::FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
106 
107 // Enable gravity
108 SET_BOOL_PROP(WaterAirBaseProblem, EnableGravity, true);
109 
110 // Use forward differences instead of central differences
111 SET_INT_PROP(WaterAirBaseProblem, NumericDifferenceMethod, +1);
112 
113 // Write newton convergence
114 SET_BOOL_PROP(WaterAirBaseProblem, NewtonWriteConvergence, false);
115 
116 // The default for the end time of the simulation (1 year)
117 SET_SCALAR_PROP(WaterAirBaseProblem, EndTime, 1.0 * 365 * 24 * 60 * 60);
118 
119 // The default for the initial time step size of the simulation
120 SET_SCALAR_PROP(WaterAirBaseProblem, InitialTimeStepSize, 250);
121 
122 // The default DGF file to load
123 SET_STRING_PROP(WaterAirBaseProblem, GridFile, "./data/waterair.dgf");
124 
125 // Use the restarted GMRES linear solver with the ILU-2 preconditioner from dune-istl
126 SET_TAG_PROP(WaterAirBaseProblem, LinearSolverSplice, ParallelIstlLinearSolver);
127 SET_TYPE_PROP(WaterAirBaseProblem, LinearSolverWrapper,
129 SET_TYPE_PROP(WaterAirBaseProblem, PreconditionerWrapper,
130  Ewoms::Linear::PreconditionerWrapperILUn<TypeTag>);
131 SET_INT_PROP(WaterAirBaseProblem, PreconditionerOrder, 2);
132 } // namespace Properties
133 } // namespace Ewoms
134 
135 namespace Ewoms {
164 template <class TypeTag >
165 class WaterAirProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
166 {
167  typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
168 
169  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
170  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
171 
172  // copy some indices for convenience
173  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
174  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
175  enum {
176  numPhases = FluidSystem::numPhases,
177 
178  // energy related indices
179  temperatureIdx = Indices::temperatureIdx,
180  energyEqIdx = Indices::energyEqIdx,
181 
182  // component indices
183  H2OIdx = FluidSystem::H2OIdx,
184  AirIdx = FluidSystem::AirIdx,
185 
186  // phase indices
187  liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
188  gasPhaseIdx = FluidSystem::gasPhaseIdx,
189 
190  // equation indices
191  conti0EqIdx = Indices::conti0EqIdx,
192 
193  // Grid and world dimension
194  dim = GridView::dimension,
195  dimWorld = GridView::dimensionworld
196  };
197 
198  static const bool enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy);
199 
200  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
201  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
202  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
203  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
204  typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;
205  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
206  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
207  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
208  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
209  typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLaw) HeatConductionLaw;
210  typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLawParams) HeatConductionLawParams;
211 
212  typedef typename GridView::ctype CoordScalar;
213  typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
214 
215  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
216 
217 public:
222  : ParentType(simulator)
223  { }
224 
228  void finishInit()
229  {
230  ParentType::finishInit();
231 
232  maxDepth_ = 1000.0; // [m]
233  eps_ = 1e-6;
234 
235  FluidSystem::init(/*Tmin=*/275, /*Tmax=*/600, /*nT=*/100,
236  /*pmin=*/9.5e6, /*pmax=*/10.5e6, /*np=*/200);
237 
238  layerBottom_ = 22.0;
239 
240  // intrinsic permeabilities
241  fineK_ = this->toDimMatrix_(1e-13);
242  coarseK_ = this->toDimMatrix_(1e-12);
243 
244  // porosities
245  finePorosity_ = 0.3;
246  coarsePorosity_ = 0.3;
247 
248  // residual saturations
249  fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
250  fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
251  coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
252  coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
253 
254  // parameters for the Brooks-Corey law
255  fineMaterialParams_.setEntryPressure(1e4);
256  coarseMaterialParams_.setEntryPressure(1e4);
257  fineMaterialParams_.setLambda(2.0);
258  coarseMaterialParams_.setLambda(2.0);
259 
260  fineMaterialParams_.finalize();
261  coarseMaterialParams_.finalize();
262 
263  // parameters for the somerton law of heat conduction
264  computeHeatCondParams_(fineHeatCondParams_, finePorosity_);
265  computeHeatCondParams_(coarseHeatCondParams_, coarsePorosity_);
266  }
267 
271 
276  std::string name() const
277  {
278  std::ostringstream oss;
279  oss << "waterair_" << Model::name();
280  if (GET_PROP_VALUE(TypeTag, EnableEnergy))
281  oss << "_ni";
282 
283  return oss.str();
284  }
285 
289  void endTimeStep()
290  {
291 #ifndef NDEBUG
292  // checkConservativeness() does not include the effect of constraints, so we
293  // disable it for this problem...
294  //this->model().checkConservativeness();
295 
296  // Calculate 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: " << storage << std::endl << std::flush;
303  }
304 #endif // NDEBUG
305  }
306 
313  template <class Context>
314  const DimMatrix& intrinsicPermeability(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
315  {
316  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
317  if (isFineMaterial_(pos))
318  return fineK_;
319  return coarseK_;
320  }
321 
325  template <class Context>
326  Scalar porosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
327  {
328  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
329  if (isFineMaterial_(pos))
330  return finePorosity_;
331  else
332  return coarsePorosity_;
333  }
334 
338  template <class Context>
339  const MaterialLawParams& materialLawParams(const Context& context,
340  unsigned spaceIdx,
341  unsigned timeIdx) const
342  {
343  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
344  if (isFineMaterial_(pos))
345  return fineMaterialParams_;
346  else
347  return coarseMaterialParams_;
348  }
349 
355  template <class Context>
356  Scalar heatCapacitySolid(const Context& context OPM_UNUSED,
357  unsigned spaceIdx OPM_UNUSED,
358  unsigned timeIdx OPM_UNUSED) const
359  {
360  return
361  790 // specific heat capacity of granite [J / (kg K)]
362  * 2700; // density of granite [kg/m^3]
363  }
364 
368  template <class Context>
369  const HeatConductionLawParams&
370  heatConductionParams(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
371  {
372  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
373  if (isFineMaterial_(pos))
374  return fineHeatCondParams_;
375  return coarseHeatCondParams_;
376  }
377 
379 
383 
393  template <class Context>
394  void boundary(BoundaryRateVector& values,
395  const Context& context,
396  unsigned spaceIdx, unsigned timeIdx) const
397  {
398  const auto& pos = context.cvCenter(spaceIdx, timeIdx);
399  assert(onLeftBoundary_(pos) ||
400  onLowerBoundary_(pos) ||
401  onRightBoundary_(pos) ||
402  onUpperBoundary_(pos));
403 
404  if (onInlet_(pos)) {
405  RateVector massRate(0.0);
406  massRate[conti0EqIdx + AirIdx] = -1e-3; // [kg/(m^2 s)]
407 
408  // impose an forced inflow boundary condition on the inlet
409  values.setMassRate(massRate);
410 
411  if (enableEnergy) {
412  Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
413  initialFluidState_(fs, context, spaceIdx, timeIdx);
414 
415  Scalar hl = fs.enthalpy(liquidPhaseIdx);
416  Scalar hg = fs.enthalpy(gasPhaseIdx);
417  values.setEnthalpyRate(values[conti0EqIdx + AirIdx] * hg +
418  values[conti0EqIdx + H2OIdx] * hl);
419  }
420  }
421  else if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
422  Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
423  initialFluidState_(fs, context, spaceIdx, timeIdx);
424 
425  // impose an freeflow boundary condition
426  values.setFreeFlow(context, spaceIdx, timeIdx, fs);
427  }
428  else
429  // no flow on top and bottom
430  values.setNoFlow();
431  }
432 
434 
438 
446  template <class Context>
447  void initial(PrimaryVariables& values,
448  const Context& context,
449  unsigned spaceIdx,
450  unsigned timeIdx) const
451  {
452  Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
453  initialFluidState_(fs, context, spaceIdx, timeIdx);
454 
455  const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
456  values.assignMassConservative(fs, matParams, /*inEquilibrium=*/true);
457  }
458 
465  template <class Context>
466  void source(RateVector& rate,
467  const Context& context OPM_UNUSED,
468  unsigned spaceIdx OPM_UNUSED,
469  unsigned timeIdx OPM_UNUSED) const
470  { rate = 0; }
471 
473 
474 private:
475  bool onLeftBoundary_(const GlobalPosition& pos) const
476  { return pos[0] < eps_; }
477 
478  bool onRightBoundary_(const GlobalPosition& pos) const
479  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
480 
481  bool onLowerBoundary_(const GlobalPosition& pos) const
482  { return pos[1] < eps_; }
483 
484  bool onUpperBoundary_(const GlobalPosition& pos) const
485  { return pos[1] > this->boundingBoxMax()[1] - eps_; }
486 
487  bool onInlet_(const GlobalPosition& pos) const
488  { return onLowerBoundary_(pos) && (15.0 < pos[0]) && (pos[0] < 25.0); }
489 
490  bool inHighTemperatureRegion_(const GlobalPosition& pos) const
491  { return (20 < pos[0]) && (pos[0] < 30) && (pos[1] < 30); }
492 
493  template <class Context, class FluidState>
494  void initialFluidState_(FluidState& fs,
495  const Context& context,
496  unsigned spaceIdx,
497  unsigned timeIdx) const
498  {
499  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
500 
501  Scalar densityW = 1000.0;
502  fs.setPressure(liquidPhaseIdx, 1e5 + (maxDepth_ - pos[1])*densityW*9.81);
503  fs.setSaturation(liquidPhaseIdx, 1.0);
504  fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
505  fs.setMoleFraction(liquidPhaseIdx, AirIdx, 0.0);
506 
507  if (inHighTemperatureRegion_(pos))
508  fs.setTemperature(380);
509  else
510  fs.setTemperature(283.0 + (maxDepth_ - pos[1])*0.03);
511 
512  // set the gas saturation and pressure
513  fs.setSaturation(gasPhaseIdx, 0);
514  Scalar pc[numPhases];
515  const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
516  MaterialLaw::capillaryPressures(pc, matParams, fs);
517  fs.setPressure(gasPhaseIdx, fs.pressure(liquidPhaseIdx) + (pc[gasPhaseIdx] - pc[liquidPhaseIdx]));
518 
519  typename FluidSystem::template ParameterCache<Scalar> paramCache;
520  typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> CFRP;
521  CFRP::solve(fs, paramCache, liquidPhaseIdx, /*setViscosity=*/false, /*setEnthalpy=*/true);
522  }
523 
524  void computeHeatCondParams_(HeatConductionLawParams& params, Scalar poro)
525  {
526  Scalar lambdaGranite = 2.8; // [W / (K m)]
527 
528  // create a Fluid state which has all phases present
529  Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
530  fs.setTemperature(293.15);
531  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
532  fs.setPressure(phaseIdx, 1.0135e5);
533  }
534 
535  typename FluidSystem::template ParameterCache<Scalar> paramCache;
536  paramCache.updateAll(fs);
537  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
538  Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
539  fs.setDensity(phaseIdx, rho);
540  }
541 
542  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
543  Scalar lambdaSaturated;
544  if (FluidSystem::isLiquid(phaseIdx)) {
545  Scalar lambdaFluid =
546  FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
547  lambdaSaturated = std::pow(lambdaGranite, (1-poro)) + std::pow(lambdaFluid, poro);
548  }
549  else
550  lambdaSaturated = std::pow(lambdaGranite, (1-poro));
551 
552  params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
553  if (!FluidSystem::isLiquid(phaseIdx))
554  params.setVacuumLambda(lambdaSaturated);
555  }
556  }
557 
558  bool isFineMaterial_(const GlobalPosition& pos) const
559  { return pos[dim-1] > layerBottom_; }
560 
561  DimMatrix fineK_;
562  DimMatrix coarseK_;
563  Scalar layerBottom_;
564 
565  Scalar finePorosity_;
566  Scalar coarsePorosity_;
567 
568  MaterialLawParams fineMaterialParams_;
569  MaterialLawParams coarseMaterialParams_;
570 
571  HeatConductionLawParams fineHeatCondParams_;
572  HeatConductionLawParams coarseHeatCondParams_;
573 
574  Scalar maxDepth_;
575  Scalar eps_;
576 };
577 } // namespace Ewoms
578 
579 #endif
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: waterairproblem.hh:228
#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: waterairproblem.hh:394
Solver wrapper for the restarted GMRES solver of dune-istl.
Definition: istlsolverwrappers.hh:125
Scalar heatCapacitySolid(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: waterairproblem.hh:356
Definition: baseauxiliarymodule.hh:37
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:339
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:326
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
const HeatConductionLawParams & heatConductionParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:370
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
#define SET_TAG_PROP(EffTypeTagName, PropTagName, ValueTypeTagName)
Define a property containing a type tag.
Definition: propertysystem.hh:436
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
Provides all unmodified linear solvers from dune-istl.
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: waterairproblem.hh:447
#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: waterairproblem.hh:276
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:314
WaterAirProblem(Simulator &simulator)
Definition: waterairproblem.hh:221
Declares the properties required for the compositional multi-phase primary variable switching model...
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
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: waterairproblem.hh:466
Non-isothermal gas injection problem where a air is injected into a fully water saturated medium...
Definition: waterairproblem.hh:57
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394
void endTimeStep()
Called by the simulator after each time integration.
Definition: waterairproblem.hh:289