My Project
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 
31 #include <opm/models/pvs/pvsproperties.hh>
32 #include <opm/simulators/linalg/parallelistlbackend.hh>
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/thermal/ConstantSolidHeatCapLaw.hpp>
42 #include <opm/material/thermal/SomertonThermalConductionLaw.hpp>
43 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.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 Opm {
56 template <class TypeTag>
57 class WaterAirProblem;
58 }
59 
60 namespace Opm::Properties {
61 
62 namespace TTag {
64 }
65 
66 // Set the grid type
67 template<class TypeTag>
68 struct Grid<TypeTag, TTag::WaterAirBaseProblem> { using type = Dune::YaspGrid<2>; };
69 
70 // Set the problem property
71 template<class TypeTag>
72 struct Problem<TypeTag, TTag::WaterAirBaseProblem> { using type = Opm::WaterAirProblem<TypeTag>; };
73 
74 // Set the material Law
75 template<class TypeTag>
76 struct MaterialLaw<TypeTag, TTag::WaterAirBaseProblem>
77 {
78 private:
79  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
80  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
81  using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
82  /*wettingPhaseIdx=*/FluidSystem::liquidPhaseIdx,
83  /*nonWettingPhaseIdx=*/FluidSystem::gasPhaseIdx>;
84 
85  // define the material law which is parameterized by effective
86  // saturations
87  using EffMaterialLaw = Opm::RegularizedBrooksCorey<Traits>;
88 
89 public:
90  // define the material law parameterized by absolute saturations
91  // which uses the two-phase API
92  using type = Opm::EffToAbsLaw<EffMaterialLaw>;
93 };
94 
95 // Set the thermal conduction law
96 template<class TypeTag>
97 struct ThermalConductionLaw<TypeTag, TTag::WaterAirBaseProblem>
98 {
99 private:
100  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
101  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
102 
103 public:
104  // define the material law parameterized by absolute saturations
105  using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
106 };
107 
108 // set the energy storage law for the solid phase
109 template<class TypeTag>
110 struct SolidEnergyLaw<TypeTag, TTag::WaterAirBaseProblem>
111 { using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
112 
113 // Set the fluid system. in this case, we use the one which describes
114 // air and water
115 template<class TypeTag>
116 struct FluidSystem<TypeTag, TTag::WaterAirBaseProblem>
117 { using type = Opm::H2OAirFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
118 
119 // Enable gravity
120 template<class TypeTag>
121 struct EnableGravity<TypeTag, TTag::WaterAirBaseProblem> { static constexpr bool value = true; };
122 
123 // Use forward differences instead of central differences
124 template<class TypeTag>
125 struct NumericDifferenceMethod<TypeTag, TTag::WaterAirBaseProblem> { static constexpr int value = +1; };
126 
127 // Write newton convergence
128 template<class TypeTag>
129 struct NewtonWriteConvergence<TypeTag, TTag::WaterAirBaseProblem> { static constexpr bool value = false; };
130 
131 // The default for the end time of the simulation (1 year)
132 template<class TypeTag>
133 struct EndTime<TypeTag, TTag::WaterAirBaseProblem>
134 {
135  using type = GetPropType<TypeTag, Scalar>;
136  static constexpr type value = 1.0 * 365 * 24 * 60 * 60;
137 };
138 
139 // The default for the initial time step size of the simulation
140 template<class TypeTag>
141 struct InitialTimeStepSize<TypeTag, TTag::WaterAirBaseProblem>
142 {
143  using type = GetPropType<TypeTag, Scalar>;
144  static constexpr type value = 250;
145 };
146 
147 // The default DGF file to load
148 template<class TypeTag>
149 struct GridFile<TypeTag, TTag::WaterAirBaseProblem> { static constexpr auto value = "./data/waterair.dgf"; };
150 
151 // Use the restarted GMRES linear solver with the ILU-2 preconditioner from dune-istl
152 template<class TypeTag>
153 struct LinearSolverSplice<TypeTag, TTag::WaterAirBaseProblem>
154 { using type = TTag::ParallelIstlLinearSolver; };
155 
156 template<class TypeTag>
157 struct LinearSolverWrapper<TypeTag, TTag::WaterAirBaseProblem>
158 { using type = Opm::Linear::SolverWrapperRestartedGMRes<TypeTag>; };
159 
160 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2,7)
161 template<class TypeTag>
162 struct PreconditionerWrapper<TypeTag, TTag::WaterAirBaseProblem>
163 { using type = Opm::Linear::PreconditionerWrapperILU<TypeTag>; };
164 #else
165 template<class TypeTag>
166 struct PreconditionerWrapper<TypeTag, TTag::WaterAirBaseProblem>
167 { using type = Opm::Linear::PreconditionerWrapperILUn<TypeTag>; };
168 #endif
169 template<class TypeTag>
170 struct PreconditionerOrder<TypeTag, TTag::WaterAirBaseProblem> { static constexpr int value = 2; };
171 
172 } // namespace Opm::Properties
173 
174 namespace Opm {
203 template <class TypeTag >
204 class WaterAirProblem : public GetPropType<TypeTag, Properties::BaseProblem>
205 {
206  using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
207 
208  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
209  using GridView = GetPropType<TypeTag, Properties::GridView>;
210 
211  // copy some indices for convenience
212  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
213  using Indices = GetPropType<TypeTag, Properties::Indices>;
214  enum {
215  numPhases = FluidSystem::numPhases,
216 
217  // energy related indices
218  temperatureIdx = Indices::temperatureIdx,
219  energyEqIdx = Indices::energyEqIdx,
220 
221  // component indices
222  H2OIdx = FluidSystem::H2OIdx,
223  AirIdx = FluidSystem::AirIdx,
224 
225  // phase indices
226  liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
227  gasPhaseIdx = FluidSystem::gasPhaseIdx,
228 
229  // equation indices
230  conti0EqIdx = Indices::conti0EqIdx,
231 
232  // Grid and world dimension
233  dim = GridView::dimension,
234  dimWorld = GridView::dimensionworld
235  };
236 
237  static const bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
238 
239  using EqVector = GetPropType<TypeTag, Properties::EqVector>;
240  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
241  using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
242  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
243  using Constraints = GetPropType<TypeTag, Properties::Constraints>;
244  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
245  using Model = GetPropType<TypeTag, Properties::Model>;
246  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
247  using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
248  using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
249  using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
250 
251  using CoordScalar = typename GridView::ctype;
252  using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
253 
254  using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
255 
256 public:
260  WaterAirProblem(Simulator& simulator)
261  : ParentType(simulator)
262  { }
263 
267  void finishInit()
268  {
269  ParentType::finishInit();
270 
271  maxDepth_ = 1000.0; // [m]
272  eps_ = 1e-6;
273 
274  FluidSystem::init(/*Tmin=*/275, /*Tmax=*/600, /*nT=*/100,
275  /*pmin=*/9.5e6, /*pmax=*/10.5e6, /*np=*/200);
276 
277  layerBottom_ = 22.0;
278 
279  // intrinsic permeabilities
280  fineK_ = this->toDimMatrix_(1e-13);
281  coarseK_ = this->toDimMatrix_(1e-12);
282 
283  // porosities
284  finePorosity_ = 0.3;
285  coarsePorosity_ = 0.3;
286 
287  // residual saturations
288  fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
289  fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
290  coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
291  coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
292 
293  // parameters for the Brooks-Corey law
294  fineMaterialParams_.setEntryPressure(1e4);
295  coarseMaterialParams_.setEntryPressure(1e4);
296  fineMaterialParams_.setLambda(2.0);
297  coarseMaterialParams_.setLambda(2.0);
298 
299  fineMaterialParams_.finalize();
300  coarseMaterialParams_.finalize();
301 
302  // parameters for the somerton law of thermal conduction
303  computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
304  computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
305 
306  // assume constant volumetric heat capacity and granite
307  solidEnergyLawParams_.setSolidHeatCapacity(790.0 // specific heat capacity of granite [J / (kg K)]
308  * 2700.0); // density of granite [kg/m^3]
309  solidEnergyLawParams_.finalize();
310  }
311 
316 
320  std::string name() const
321  {
322  std::ostringstream oss;
323  oss << "waterair_" << Model::name();
324  if (getPropValue<TypeTag, Properties::EnableEnergy>())
325  oss << "_ni";
326 
327  return oss.str();
328  }
329 
333  void endTimeStep()
334  {
335 #ifndef NDEBUG
336  // checkConservativeness() does not include the effect of constraints, so we
337  // disable it for this problem...
338  //this->model().checkConservativeness();
339 
340  // Calculate storage terms
341  EqVector storage;
342  this->model().globalStorage(storage);
343 
344  // Write mass balance information for rank 0
345  if (this->gridView().comm().rank() == 0) {
346  std::cout << "Storage: " << storage << std::endl << std::flush;
347  }
348 #endif // NDEBUG
349  }
350 
357  template <class Context>
358  const DimMatrix& intrinsicPermeability(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
359  {
360  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
361  if (isFineMaterial_(pos))
362  return fineK_;
363  return coarseK_;
364  }
365 
369  template <class Context>
370  Scalar porosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
371  {
372  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
373  if (isFineMaterial_(pos))
374  return finePorosity_;
375  else
376  return coarsePorosity_;
377  }
378 
382  template <class Context>
383  const MaterialLawParams& materialLawParams(const Context& context,
384  unsigned spaceIdx,
385  unsigned timeIdx) const
386  {
387  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
388  if (isFineMaterial_(pos))
389  return fineMaterialParams_;
390  else
391  return coarseMaterialParams_;
392  }
393 
399  template <class Context>
400  const SolidEnergyLawParams&
401  solidEnergyLawParams(const Context& /*context*/,
402  unsigned /*spaceIdx*/,
403  unsigned /*timeIdx*/) const
404  { return solidEnergyLawParams_; }
405 
409  template <class Context>
410  const ThermalConductionLawParams&
411  thermalConductionLawParams(const Context& context,
412  unsigned spaceIdx,
413  unsigned timeIdx) const
414  {
415  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
416  if (isFineMaterial_(pos))
417  return fineThermalCondParams_;
418  return coarseThermalCondParams_;
419  }
420 
422 
427 
436  template <class Context>
437  void boundary(BoundaryRateVector& values,
438  const Context& context,
439  unsigned spaceIdx, unsigned timeIdx) const
440  {
441  const auto& pos = context.cvCenter(spaceIdx, timeIdx);
442  assert(onLeftBoundary_(pos) ||
443  onLowerBoundary_(pos) ||
444  onRightBoundary_(pos) ||
445  onUpperBoundary_(pos));
446 
447  if (onInlet_(pos)) {
448  RateVector massRate(0.0);
449  massRate[conti0EqIdx + AirIdx] = -1e-3; // [kg/(m^2 s)]
450 
451  // impose an forced inflow boundary condition on the inlet
452  values.setMassRate(massRate);
453 
454  if (enableEnergy) {
455  Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
456  initialFluidState_(fs, context, spaceIdx, timeIdx);
457 
458  Scalar hl = fs.enthalpy(liquidPhaseIdx);
459  Scalar hg = fs.enthalpy(gasPhaseIdx);
460  values.setEnthalpyRate(values[conti0EqIdx + AirIdx] * hg +
461  values[conti0EqIdx + H2OIdx] * hl);
462  }
463  }
464  else if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
465  Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
466  initialFluidState_(fs, context, spaceIdx, timeIdx);
467 
468  // impose an freeflow boundary condition
469  values.setFreeFlow(context, spaceIdx, timeIdx, fs);
470  }
471  else
472  // no flow on top and bottom
473  values.setNoFlow();
474  }
475 
477 
482 
489  template <class Context>
490  void initial(PrimaryVariables& values,
491  const Context& context,
492  unsigned spaceIdx,
493  unsigned timeIdx) const
494  {
495  Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
496  initialFluidState_(fs, context, spaceIdx, timeIdx);
497 
498  const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
499  values.assignMassConservative(fs, matParams, /*inEquilibrium=*/true);
500  }
501 
508  template <class Context>
509  void source(RateVector& rate,
510  const Context& /*context*/,
511  unsigned /*spaceIdx*/,
512  unsigned /*timeIdx*/) const
513  { rate = 0; }
514 
516 
517 private:
518  bool onLeftBoundary_(const GlobalPosition& pos) const
519  { return pos[0] < eps_; }
520 
521  bool onRightBoundary_(const GlobalPosition& pos) const
522  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
523 
524  bool onLowerBoundary_(const GlobalPosition& pos) const
525  { return pos[1] < eps_; }
526 
527  bool onUpperBoundary_(const GlobalPosition& pos) const
528  { return pos[1] > this->boundingBoxMax()[1] - eps_; }
529 
530  bool onInlet_(const GlobalPosition& pos) const
531  { return onLowerBoundary_(pos) && (15.0 < pos[0]) && (pos[0] < 25.0); }
532 
533  bool inHighTemperatureRegion_(const GlobalPosition& pos) const
534  { return (20 < pos[0]) && (pos[0] < 30) && (pos[1] < 30); }
535 
536  template <class Context, class FluidState>
537  void initialFluidState_(FluidState& fs,
538  const Context& context,
539  unsigned spaceIdx,
540  unsigned timeIdx) const
541  {
542  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
543 
544  Scalar densityW = 1000.0;
545  fs.setPressure(liquidPhaseIdx, 1e5 + (maxDepth_ - pos[1])*densityW*9.81);
546  fs.setSaturation(liquidPhaseIdx, 1.0);
547  fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
548  fs.setMoleFraction(liquidPhaseIdx, AirIdx, 0.0);
549 
550  if (inHighTemperatureRegion_(pos))
551  fs.setTemperature(380);
552  else
553  fs.setTemperature(283.0 + (maxDepth_ - pos[1])*0.03);
554 
555  // set the gas saturation and pressure
556  fs.setSaturation(gasPhaseIdx, 0);
557  Scalar pc[numPhases];
558  const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
559  MaterialLaw::capillaryPressures(pc, matParams, fs);
560  fs.setPressure(gasPhaseIdx, fs.pressure(liquidPhaseIdx) + (pc[gasPhaseIdx] - pc[liquidPhaseIdx]));
561 
562  typename FluidSystem::template ParameterCache<Scalar> paramCache;
563  using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
564  CFRP::solve(fs, paramCache, liquidPhaseIdx, /*setViscosity=*/true, /*setEnthalpy=*/true);
565  }
566 
567  void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
568  {
569  Scalar lambdaGranite = 2.8; // [W / (K m)]
570 
571  // create a Fluid state which has all phases present
572  Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
573  fs.setTemperature(293.15);
574  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
575  fs.setPressure(phaseIdx, 1.0135e5);
576  }
577 
578  typename FluidSystem::template ParameterCache<Scalar> paramCache;
579  paramCache.updateAll(fs);
580  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
581  Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
582  fs.setDensity(phaseIdx, rho);
583  }
584 
585  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
586  Scalar lambdaSaturated;
587  if (FluidSystem::isLiquid(phaseIdx)) {
588  Scalar lambdaFluid =
589  FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
590  lambdaSaturated = std::pow(lambdaGranite, (1-poro)) + std::pow(lambdaFluid, poro);
591  }
592  else
593  lambdaSaturated = std::pow(lambdaGranite, (1-poro));
594 
595  params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
596  if (!FluidSystem::isLiquid(phaseIdx))
597  params.setVacuumLambda(lambdaSaturated);
598  }
599  }
600 
601  bool isFineMaterial_(const GlobalPosition& pos) const
602  { return pos[dim-1] > layerBottom_; }
603 
604  DimMatrix fineK_;
605  DimMatrix coarseK_;
606  Scalar layerBottom_;
607 
608  Scalar finePorosity_;
609  Scalar coarsePorosity_;
610 
611  MaterialLawParams fineMaterialParams_;
612  MaterialLawParams coarseMaterialParams_;
613 
614  ThermalConductionLawParams fineThermalCondParams_;
615  ThermalConductionLawParams coarseThermalCondParams_;
616  SolidEnergyLawParams solidEnergyLawParams_;
617 
618  Scalar maxDepth_;
619  Scalar eps_;
620 };
621 } // namespace Opm
622 
623 #endif
Non-isothermal gas injection problem where a air is injected into a fully water saturated medium.
Definition: waterairproblem.hh:205
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:358
const SolidEnergyLawParams & solidEnergyLawParams(const Context &, unsigned, unsigned) const
Return the parameters for the energy storage law of the rock.
Definition: waterairproblem.hh:401
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:411
void finishInit()
Definition: waterairproblem.hh:267
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:490
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:437
std::string name() const
Definition: waterairproblem.hh:320
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:370
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: waterairproblem.hh:509
void endTimeStep()
Definition: waterairproblem.hh:333
WaterAirProblem(Simulator &simulator)
Definition: waterairproblem.hh:260
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:383
Definition: waterairproblem.hh:63