My Project
Loading...
Searching...
No Matches
powerinjectionproblem.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_POWER_INJECTION_PROBLEM_HH
29#define EWOMS_POWER_INJECTION_PROBLEM_HH
30
31#include <opm/models/immiscible/immisciblemodel.hh>
32#include <opm/models/io/cubegridvanguard.hh>
33
34#include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
35#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
36#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
37#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
38#include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp>
39#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
40#include <opm/material/components/SimpleH2O.hpp>
41#include <opm/material/components/Air.hpp>
42
43#include <dune/grid/yaspgrid.hh>
44
45#include <dune/common/version.hh>
46#include <dune/common/fvector.hh>
47#include <dune/common/fmatrix.hh>
48
49#include <sstream>
50#include <string>
51#include <type_traits>
52#include <iostream>
53
54namespace Opm {
55template <class TypeTag>
56class PowerInjectionProblem;
57}
58
59namespace Opm::Properties {
60
61namespace TTag {
63}
64
65// Set the grid implementation to be used
66template<class TypeTag>
67struct Grid<TypeTag, TTag::PowerInjectionBaseProblem> { using type = Dune::YaspGrid</*dim=*/1>; };
68
69// set the Vanguard property
70template<class TypeTag>
71struct Vanguard<TypeTag, TTag::PowerInjectionBaseProblem> { using type = Opm::CubeGridVanguard<TypeTag>; };
72
73// Set the problem property
74template<class TypeTag>
75struct Problem<TypeTag, TTag::PowerInjectionBaseProblem> { using type = Opm::PowerInjectionProblem<TypeTag>; };
76
77// Set the wetting phase
78template<class TypeTag>
79struct WettingPhase<TypeTag, TTag::PowerInjectionBaseProblem>
80{
81private:
82 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
83
84public:
85 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
86};
87
88// Set the non-wetting phase
89template<class TypeTag>
90struct NonwettingPhase<TypeTag, TTag::PowerInjectionBaseProblem>
91{
92private:
93 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
94
95public:
96 using type = Opm::GasPhase<Scalar, Opm::Air<Scalar> >;
97};
98
99// Set the material Law
100template<class TypeTag>
101struct MaterialLaw<TypeTag, TTag::PowerInjectionBaseProblem>
102{
103private:
104 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
105 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
106 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
107
108 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
109 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
110 /*wettingPhaseIdx=*/FluidSystem::wettingPhaseIdx,
111 /*nonWettingPhaseIdx=*/FluidSystem::nonWettingPhaseIdx>;
112
113 // define the material law which is parameterized by effective
114 // saturations
115 using EffectiveLaw = Opm::RegularizedVanGenuchten<Traits>;
116
117public:
118 // define the material law parameterized by absolute saturations
119 using type = Opm::EffToAbsLaw<EffectiveLaw>;
120};
121
122// Write out the filter velocities for this problem
123template<class TypeTag>
124struct VtkWriteFilterVelocities<TypeTag, TTag::PowerInjectionBaseProblem> { static constexpr bool value = true; };
125
126// Disable gravity
127template<class TypeTag>
128struct EnableGravity<TypeTag, TTag::PowerInjectionBaseProblem> { static constexpr bool value = false; };
129
130// define the properties specific for the power injection problem
131template<class TypeTag>
132struct DomainSizeX<TypeTag, TTag::PowerInjectionBaseProblem>
133{
134 using type = GetPropType<TypeTag, Scalar>;
135 static constexpr type value = 100.0;
136};
137template<class TypeTag>
138struct DomainSizeY<TypeTag, TTag::PowerInjectionBaseProblem>
139{
140 using type = GetPropType<TypeTag, Scalar>;
141 static constexpr type value = 1.0;
142};
143template<class TypeTag>
144struct DomainSizeZ<TypeTag, TTag::PowerInjectionBaseProblem>
145{
146 using type = GetPropType<TypeTag, Scalar>;
147 static constexpr type value = 1.0;
148};
149
150template<class TypeTag>
151struct CellsX<TypeTag, TTag::PowerInjectionBaseProblem> { static constexpr int value = 250; };
152template<class TypeTag>
153struct CellsY<TypeTag, TTag::PowerInjectionBaseProblem> { static constexpr int value = 1; };
154template<class TypeTag>
155struct CellsZ<TypeTag, TTag::PowerInjectionBaseProblem> { static constexpr int value = 1; };
156
157// The default for the end time of the simulation
158template<class TypeTag>
159struct EndTime<TypeTag, TTag::PowerInjectionBaseProblem>
160{
161 using type = GetPropType<TypeTag, Scalar>;
162 static constexpr type value = 100;
163};
164
165// The default for the initial time step size of the simulation
166template<class TypeTag>
167struct InitialTimeStepSize<TypeTag, TTag::PowerInjectionBaseProblem>
168{
169 using type = GetPropType<TypeTag, Scalar>;
170 static constexpr type value = 1e-3;
171};
172
173} // namespace Opm::Properties
174
175namespace Opm {
188template <class TypeTag>
189class PowerInjectionProblem : public GetPropType<TypeTag, Properties::BaseProblem>
190{
191 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
192
193 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
194 using GridView = GetPropType<TypeTag, Properties::GridView>;
195 using Indices = GetPropType<TypeTag, Properties::Indices>;
196 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
197 using WettingPhase = GetPropType<TypeTag, Properties::WettingPhase>;
198 using NonwettingPhase = GetPropType<TypeTag, Properties::NonwettingPhase>;
199 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
200 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
201 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
202 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
203 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
204
205 enum {
206 // number of phases
207 numPhases = FluidSystem::numPhases,
208
209 // phase indices
210 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
211 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
212
213 // equation indices
214 contiNEqIdx = Indices::conti0EqIdx + nonWettingPhaseIdx,
215
216 // Grid and world dimension
217 dim = GridView::dimension,
218 dimWorld = GridView::dimensionworld
219 };
220
221 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
222 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
223
224 using CoordScalar = typename GridView::ctype;
225 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
226
227 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
228
229public:
233 PowerInjectionProblem(Simulator& simulator)
234 : ParentType(simulator)
235 { }
236
241 {
242 ParentType::finishInit();
243
244 eps_ = 3e-6;
245 FluidSystem::init();
246
247 temperature_ = 273.15 + 26.6;
248
249 // parameters for the Van Genuchten law
250 // alpha and n
251 materialParams_.setVgAlpha(0.00045);
252 materialParams_.setVgN(7.3);
253 materialParams_.finalize();
254
255 K_ = this->toDimMatrix_(5.73e-08); // [m^2]
256
257 setupInitialFluidState_();
258 }
259
264
268 std::string name() const
269 {
270 std::ostringstream oss;
271 oss << "powerinjection_";
272 if (std::is_same<GetPropType<TypeTag, Properties::FluxModule>,
273 Opm::DarcyFluxModule<TypeTag> >::value)
274 oss << "darcy";
275 else
276 oss << "forchheimer";
277
278 if (std::is_same<GetPropType<TypeTag, Properties::LocalLinearizerSplice>,
279 Properties::TTag::AutoDiffLocalLinearizer>::value)
280 oss << "_" << "ad";
281 else
282 oss << "_" << "fd";
283
284 return oss.str();
285 }
286
291 {
292#ifndef NDEBUG
293 this->model().checkConservativeness();
294
295 // Calculate storage terms
296 EqVector storage;
297 this->model().globalStorage(storage);
298
299 // Write mass balance information for rank 0
300 if (this->gridView().comm().rank() == 0) {
301 std::cout << "Storage: " << storage << std::endl << std::flush;
302 }
303#endif // NDEBUG
304 }
306
311
315 template <class Context>
316 const DimMatrix& intrinsicPermeability(const Context& /*context*/,
317 unsigned /*spaceIdx*/,
318 unsigned /*timeIdx*/) const
319 { return K_; }
320
324 template <class Context>
325 Scalar ergunCoefficient(const Context& /*context*/,
326 unsigned /*spaceIdx*/,
327 unsigned /*timeIdx*/) const
328 { return 0.3866; }
329
333 template <class Context>
334 Scalar porosity(const Context& /*context*/,
335 unsigned /*spaceIdx*/,
336 unsigned /*timeIdx*/) const
337 { return 0.558; }
338
342 template <class Context>
343 const MaterialLawParams&
344 materialLawParams(const Context& /*context*/,
345 unsigned /*spaceIdx*/,
346 unsigned /*timeIdx*/) const
347 { return materialParams_; }
348
352 template <class Context>
353 Scalar temperature(const Context& /*context*/,
354 unsigned /*spaceIdx*/,
355 unsigned /*timeIdx*/) const
356 { return temperature_; }
357
359
364
371 template <class Context>
372 void boundary(BoundaryRateVector& values,
373 const Context& context,
374 unsigned spaceIdx,
375 unsigned timeIdx) const
376 {
377 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
378
379 if (onLeftBoundary_(pos)) {
380 RateVector massRate(0.0);
381 massRate = 0.0;
382 massRate[contiNEqIdx] = -1.00; // kg / (m^2 * s)
383
384 // impose a forced flow boundary
385 values.setMassRate(massRate);
386 }
387 else if (onRightBoundary_(pos))
388 // free flow boundary with initial condition on the right
389 values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
390 else
391 values.setNoFlow();
392 }
393
395
400
404 template <class Context>
405 void initial(PrimaryVariables& values,
406 const Context& /*context*/,
407 unsigned /*spaceIdx*/,
408 unsigned /*timeIdx*/) const
409 {
410 // assign the primary variables
411 values.assignNaive(initialFluidState_);
412 }
413
420 template <class Context>
421 void source(RateVector& rate,
422 const Context& /*context*/,
423 unsigned /*spaceIdx*/,
424 unsigned /*timeIdx*/) const
425 { rate = Scalar(0.0); }
426
428
429private:
430 bool onLeftBoundary_(const GlobalPosition& pos) const
431 { return pos[0] < this->boundingBoxMin()[0] + eps_; }
432
433 bool onRightBoundary_(const GlobalPosition& pos) const
434 { return pos[0] > this->boundingBoxMax()[0] - eps_; }
435
436 void setupInitialFluidState_()
437 {
438 initialFluidState_.setTemperature(temperature_);
439
440 Scalar Sw = 1.0;
441 initialFluidState_.setSaturation(wettingPhaseIdx, Sw);
442 initialFluidState_.setSaturation(nonWettingPhaseIdx, 1 - Sw);
443
444 Scalar p = 1e5;
445 initialFluidState_.setPressure(wettingPhaseIdx, p);
446 initialFluidState_.setPressure(nonWettingPhaseIdx, p);
447
448 typename FluidSystem::template ParameterCache<Scalar> paramCache;
449 paramCache.updateAll(initialFluidState_);
450 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
451 initialFluidState_.setDensity(phaseIdx,
452 FluidSystem::density(initialFluidState_, paramCache, phaseIdx));
453 initialFluidState_.setViscosity(phaseIdx,
454 FluidSystem::viscosity(initialFluidState_, paramCache, phaseIdx));
455 }
456 }
457
458 DimMatrix K_;
459 MaterialLawParams materialParams_;
460
461 Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;
462 Scalar temperature_;
463 Scalar eps_;
464};
465
466} // namespace Opm
467
468#endif
1D Problem with very fast injection of gas on the left.
Definition powerinjectionproblem.hh:190
void finishInit()
Definition powerinjectionproblem.hh:240
Scalar porosity(const Context &, unsigned, unsigned) const
Definition powerinjectionproblem.hh:334
Scalar temperature(const Context &, unsigned, unsigned) const
Definition powerinjectionproblem.hh:353
void endTimeStep()
Definition powerinjectionproblem.hh:290
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Definition powerinjectionproblem.hh:405
Scalar ergunCoefficient(const Context &, unsigned, unsigned) const
Definition powerinjectionproblem.hh:325
const MaterialLawParams & materialLawParams(const Context &, unsigned, unsigned) const
Definition powerinjectionproblem.hh:344
std::string name() const
Definition powerinjectionproblem.hh:268
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition powerinjectionproblem.hh:421
PowerInjectionProblem(Simulator &simulator)
Definition powerinjectionproblem.hh:233
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition powerinjectionproblem.hh:372
const DimMatrix & intrinsicPermeability(const Context &, unsigned, unsigned) const
Definition powerinjectionproblem.hh:316
Definition powerinjectionproblem.hh:62