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