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
55namespace Opm {
56template <class TypeTag>
57class WaterAirProblem;
58}
59
60namespace Opm::Properties {
61
62namespace TTag {
64}
65
66// Set the grid type
67template<class TypeTag>
68struct Grid<TypeTag, TTag::WaterAirBaseProblem> { using type = Dune::YaspGrid<2>; };
69
70// Set the problem property
71template<class TypeTag>
72struct Problem<TypeTag, TTag::WaterAirBaseProblem> { using type = Opm::WaterAirProblem<TypeTag>; };
73
74// Set the material Law
75template<class TypeTag>
76struct MaterialLaw<TypeTag, TTag::WaterAirBaseProblem>
77{
78private:
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
89public:
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
96template<class TypeTag>
97struct ThermalConductionLaw<TypeTag, TTag::WaterAirBaseProblem>
98{
99private:
100 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
101 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
102
103public:
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
109template<class TypeTag>
110struct 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
115template<class TypeTag>
116struct FluidSystem<TypeTag, TTag::WaterAirBaseProblem>
117{ using type = Opm::H2OAirFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
118
119// Enable gravity
120template<class TypeTag>
121struct EnableGravity<TypeTag, TTag::WaterAirBaseProblem> { static constexpr bool value = true; };
122
123// Use forward differences instead of central differences
124template<class TypeTag>
125struct NumericDifferenceMethod<TypeTag, TTag::WaterAirBaseProblem> { static constexpr int value = +1; };
126
127// Write newton convergence
128template<class TypeTag>
129struct NewtonWriteConvergence<TypeTag, TTag::WaterAirBaseProblem> { static constexpr bool value = false; };
130
131// The default for the end time of the simulation (1 year)
132template<class TypeTag>
133struct 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
140template<class TypeTag>
141struct 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
148template<class TypeTag>
149struct 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
152template<class TypeTag>
153struct LinearSolverSplice<TypeTag, TTag::WaterAirBaseProblem>
154{ using type = TTag::ParallelIstlLinearSolver; };
155
156template<class TypeTag>
157struct LinearSolverWrapper<TypeTag, TTag::WaterAirBaseProblem>
158{ using type = Opm::Linear::SolverWrapperRestartedGMRes<TypeTag>; };
159
160#if DUNE_VERSION_NEWER(DUNE_ISTL, 2,7)
161template<class TypeTag>
162struct PreconditionerWrapper<TypeTag, TTag::WaterAirBaseProblem>
163{ using type = Opm::Linear::PreconditionerWrapperILU<TypeTag>; };
164#else
165template<class TypeTag>
166struct PreconditionerWrapper<TypeTag, TTag::WaterAirBaseProblem>
167{ using type = Opm::Linear::PreconditionerWrapperILUn<TypeTag>; };
168#endif
169template<class TypeTag>
170struct PreconditionerOrder<TypeTag, TTag::WaterAirBaseProblem> { static constexpr int value = 2; };
171
172} // namespace Opm::Properties
173
174namespace Opm {
203template <class TypeTag >
204class 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
256public:
260 WaterAirProblem(Simulator& simulator)
261 : ParentType(simulator)
262 { }
263
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
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
517private:
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
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
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: waterairproblem.hh:411
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
const SolidEnergyLawParams & solidEnergyLawParams(const Context &, unsigned, unsigned) const
Return the parameters for the energy storage law of the rock.
Definition: waterairproblem.hh:401
Definition: waterairproblem.hh:63