My Project
Loading...
Searching...
No Matches
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
54namespace Opm {
55template <class TypeTag>
56class ObstacleProblem;
57}
58
59namespace Opm::Properties {
60
61namespace TTag {
63}
64
65// Set the grid type
66template<class TypeTag>
67struct Grid<TypeTag, TTag::ObstacleBaseProblem> { using type = Dune::YaspGrid<2>; };
68
69// Set the problem property
70template<class TypeTag>
71struct Problem<TypeTag, TTag::ObstacleBaseProblem> { using type = Opm::ObstacleProblem<TypeTag>; };
72
73// Set fluid configuration
74template<class TypeTag>
75struct FluidSystem<TypeTag, TTag::ObstacleBaseProblem>
76{ using type = Opm::H2ON2FluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
77
78// Set the material Law
79template<class TypeTag>
80struct MaterialLaw<TypeTag, TTag::ObstacleBaseProblem>
81{
82private:
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
92public:
93 using type = Opm::EffToAbsLaw<EffMaterialLaw>;
94};
95
96// Set the thermal conduction law
97template<class TypeTag>
98struct ThermalConductionLaw<TypeTag, TTag::ObstacleBaseProblem>
99{
100private:
101 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
102 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
103
104public:
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
110template<class TypeTag>
111struct SolidEnergyLaw<TypeTag, TTag::ObstacleBaseProblem>
112{ using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
113
114// Enable gravity
115template<class TypeTag>
116struct EnableGravity<TypeTag, TTag::ObstacleBaseProblem> { static constexpr bool value = true; };
117
118// The default for the end time of the simulation
119template<class TypeTag>
120struct 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
127template<class TypeTag>
128struct 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
135template<class TypeTag>
136struct GridFile<TypeTag, TTag::ObstacleBaseProblem> { static constexpr auto value = "./data/obstacle_24x16.dgf"; };
137
138} // namespace Opm::Properties
139
140namespace Opm {
167template <class TypeTag>
168class 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
201public:
205 ObstacleProblem(Simulator& simulator)
206 : ParentType(simulator)
207 { }
208
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
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
468private:
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
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
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition obstacleproblem.hh:340
const ThermalConductionLawParams & thermalConductionParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition obstacleproblem.hh:397
Scalar temperature(const Context &, unsigned, unsigned) const
Definition obstacleproblem.hh:330
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
const SolidEnergyLawParams & solidEnergyLawParams(const Context &, unsigned, unsigned) const
Return the parameters for the energy storage law of the rock.
Definition obstacleproblem.hh:387
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition obstacleproblem.hh:460
ObstacleProblem(Simulator &simulator)
Definition obstacleproblem.hh:205
Definition obstacleproblem.hh:62