28 #ifndef EWOMS_CUVETTE_PROBLEM_HH
29 #define EWOMS_CUVETTE_PROBLEM_HH
31 #include <opm/models/pvs/pvsproperties.hh>
33 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
34 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
35 #include <opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp>
36 #include <opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp>
37 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
38 #include <opm/material/thermal/ConstantSolidHeatCapLaw.hpp>
39 #include <opm/material/thermal/SomertonThermalConductionLaw.hpp>
40 #include <opm/material/constraintsolvers/MiscibleMultiPhaseComposition.hpp>
41 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
42 #include <opm/material/common/Valgrind.hpp>
43 #include <opm/material/common/Unused.hpp>
45 #include <dune/grid/yaspgrid.hh>
46 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
48 #include <dune/common/version.hh>
49 #include <dune/common/fvector.hh>
50 #include <dune/common/fmatrix.hh>
55 template <
class TypeTag>
59 namespace Opm::Properties {
68 template<
class TypeTag>
69 struct Grid<TypeTag, TTag::CuvetteBaseProblem> {
using type = Dune::YaspGrid<2>; };
72 template<
class TypeTag>
76 template<
class TypeTag>
77 struct FluidSystem<TypeTag, TTag::CuvetteBaseProblem>
78 {
using type = Opm::H2OAirMesityleneFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
81 template<
class TypeTag>
82 struct EnableGravity<TypeTag, TTag::CuvetteBaseProblem> {
static constexpr
bool value =
true; };
85 template<
class TypeTag>
86 struct MaxTimeStepSize<TypeTag, TTag::CuvetteBaseProblem>
88 using type = GetPropType<TypeTag, Scalar>;
89 static constexpr type value = 600.;
93 template<
class TypeTag>
94 struct MaterialLaw<TypeTag, TTag::CuvetteBaseProblem>
97 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
98 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
100 using Traits = Opm::ThreePhaseMaterialTraits<
102 FluidSystem::waterPhaseIdx,
103 FluidSystem::naplPhaseIdx,
104 FluidSystem::gasPhaseIdx>;
107 using type = Opm::ThreePhaseParkerVanGenuchten<Traits>;
111 template<
class TypeTag>
112 struct SolidEnergyLaw<TypeTag, TTag::CuvetteBaseProblem>
113 {
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
116 template<
class TypeTag>
117 struct ThermalConductionLaw<TypeTag, TTag::CuvetteBaseProblem>
120 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
121 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
125 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
129 template<
class TypeTag>
130 struct EndTime<TypeTag, TTag::CuvetteBaseProblem>
132 using type = GetPropType<TypeTag, Scalar>;
133 static constexpr type value = 180;
137 template<
class TypeTag>
138 struct InitialTimeStepSize<TypeTag, TTag::CuvetteBaseProblem>
140 using type = GetPropType<TypeTag, Scalar>;
141 static constexpr type value = 1;
145 template<
class TypeTag>
146 struct GridFile<TypeTag, TTag::CuvetteBaseProblem> {
static constexpr
auto value =
"./data/cuvette_11x4.dgf"; };
178 template <
class TypeTag>
181 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
183 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
184 using GridView = GetPropType<TypeTag, Properties::GridView>;
185 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
186 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
187 using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
188 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
189 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
190 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
191 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
192 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
193 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
194 using Model = GetPropType<TypeTag, Properties::Model>;
195 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
198 using Indices = GetPropType<TypeTag, Properties::Indices>;
199 enum { numPhases = FluidSystem::numPhases };
200 enum { numComponents = FluidSystem::numComponents };
201 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
202 enum { naplPhaseIdx = FluidSystem::naplPhaseIdx };
203 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
204 enum { H2OIdx = FluidSystem::H2OIdx };
205 enum { airIdx = FluidSystem::airIdx };
206 enum { NAPLIdx = FluidSystem::NAPLIdx };
207 enum { conti0EqIdx = Indices::conti0EqIdx };
210 enum { dimWorld = GridView::dimensionworld };
212 using CoordScalar =
typename GridView::ctype;
213 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
214 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
221 : ParentType(simulator)
230 ParentType::finishInit();
232 if (Opm::Valgrind::IsRunning())
233 FluidSystem::init(283.15, 500.0, 20,
236 FluidSystem::init(283.15, 500.0, 200,
240 fineK_ = this->toDimMatrix_(6.28e-12);
241 coarseK_ = this->toDimMatrix_(9.14e-10);
244 finePorosity_ = 0.42;
245 coarsePorosity_ = 0.42;
250 fineMaterialParams_.setVgAlpha(0.0005);
251 coarseMaterialParams_.setVgAlpha(0.005);
252 fineMaterialParams_.setVgN(4.0);
253 coarseMaterialParams_.setVgN(4.0);
255 coarseMaterialParams_.setkrRegardsSnr(
true);
256 fineMaterialParams_.setkrRegardsSnr(
true);
259 fineMaterialParams_.setSwr(0.1201);
260 fineMaterialParams_.setSwrx(0.1201);
261 fineMaterialParams_.setSnr(0.0701);
262 fineMaterialParams_.setSgr(0.0101);
263 coarseMaterialParams_.setSwr(0.1201);
264 coarseMaterialParams_.setSwrx(0.1201);
265 coarseMaterialParams_.setSnr(0.0701);
266 coarseMaterialParams_.setSgr(0.0101);
269 fineMaterialParams_.setPcMinSat(gasPhaseIdx, 0);
270 fineMaterialParams_.setPcMaxSat(gasPhaseIdx, 0);
271 fineMaterialParams_.setPcMinSat(naplPhaseIdx, 0);
272 fineMaterialParams_.setPcMaxSat(naplPhaseIdx, -1000);
273 fineMaterialParams_.setPcMinSat(waterPhaseIdx, 0);
274 fineMaterialParams_.setPcMaxSat(waterPhaseIdx, -10000);
276 coarseMaterialParams_.setPcMinSat(gasPhaseIdx, 0);
277 coarseMaterialParams_.setPcMaxSat(gasPhaseIdx, 0);
278 coarseMaterialParams_.setPcMinSat(naplPhaseIdx, 0);
279 coarseMaterialParams_.setPcMaxSat(naplPhaseIdx, -100);
280 coarseMaterialParams_.setPcMinSat(waterPhaseIdx, 0);
281 coarseMaterialParams_.setPcMaxSat(waterPhaseIdx, -1000);
284 fineMaterialParams_.setResidSat(waterPhaseIdx, 0.1201);
285 fineMaterialParams_.setResidSat(naplPhaseIdx, 0.0701);
286 fineMaterialParams_.setResidSat(gasPhaseIdx, 0.0101);
288 coarseMaterialParams_.setResidSat(waterPhaseIdx, 0.1201);
289 coarseMaterialParams_.setResidSat(naplPhaseIdx, 0.0701);
290 coarseMaterialParams_.setResidSat(gasPhaseIdx, 0.0101);
293 fineMaterialParams_.finalize();
294 coarseMaterialParams_.finalize();
297 computeThermalCondParams_(thermalCondParams_, finePorosity_);
300 solidEnergyLawParams_.setSolidHeatCapacity(790.0
302 solidEnergyLawParams_.finalize();
304 initInjectFluidState_();
324 {
return std::string(
"cuvette_") + Model::name(); }
332 this->model().checkConservativeness();
336 this->model().globalStorage(storage);
339 if (this->gridView().comm().rank() == 0) {
340 std::cout <<
"Storage: " << storage << std::endl << std::flush;
355 template <
class Context>
357 unsigned spaceIdx OPM_UNUSED,
358 unsigned timeIdx OPM_UNUSED)
const
364 template <
class Context>
366 unsigned timeIdx)
const
368 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
369 if (isFineMaterial_(pos))
377 template <
class Context>
378 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
380 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
381 if (isFineMaterial_(pos))
382 return finePorosity_;
384 return coarsePorosity_;
390 template <
class Context>
392 unsigned spaceIdx,
unsigned timeIdx)
const
394 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
395 if (isFineMaterial_(pos))
396 return fineMaterialParams_;
398 return coarseMaterialParams_;
404 template <
class Context>
405 const ThermalConductionLawParams &
407 unsigned spaceIdx OPM_UNUSED,
408 unsigned timeIdx OPM_UNUSED)
const
409 {
return thermalCondParams_; }
421 template <
class Context>
422 void boundary(BoundaryRateVector& values,
const Context& context,
423 unsigned spaceIdx,
unsigned timeIdx)
const
425 const auto& pos = context.pos(spaceIdx, timeIdx);
427 if (onRightBoundary_(pos)) {
428 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
430 initialFluidState_(fs, context, spaceIdx, timeIdx);
432 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
435 else if (onLeftBoundary_(pos)) {
437 RateVector molarRate;
441 Scalar molarInjectionRate = 0.3435;
442 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
443 molarRate[conti0EqIdx + compIdx] =
445 * injectFluidState_.moleFraction(gasPhaseIdx, compIdx);
448 Scalar massInjectionRate =
450 * injectFluidState_.averageMolarMass(gasPhaseIdx);
453 values.setMolarRate(molarRate);
454 values.setEnthalpyRate(-injectFluidState_.enthalpy(gasPhaseIdx) * massInjectionRate);
470 template <
class Context>
471 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
472 unsigned timeIdx)
const
474 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
476 initialFluidState_(fs, context, spaceIdx, timeIdx);
479 values.assignMassConservative(fs, matParams,
false);
488 template <
class Context>
490 const Context& context OPM_UNUSED,
491 unsigned spaceIdx OPM_UNUSED,
492 unsigned timeIdx OPM_UNUSED)
const
493 { rate = Scalar(0.0); }
498 bool onLeftBoundary_(
const GlobalPosition& pos)
const
499 {
return pos[0] < eps_; }
501 bool onRightBoundary_(
const GlobalPosition& pos)
const
502 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
504 bool onLowerBoundary_(
const GlobalPosition& pos)
const
505 {
return pos[1] < eps_; }
507 bool onUpperBoundary_(
const GlobalPosition& pos)
const
508 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
510 bool isContaminated_(
const GlobalPosition& pos)
const
512 return (0.20 <= pos[0]) && (pos[0] <= 0.80) && (0.4 <= pos[1])
516 bool isFineMaterial_(
const GlobalPosition& pos)
const
518 if (0.13 <= pos[0] && 1.20 >= pos[0] && 0.32 <= pos[1] && pos[1] <= 0.57)
520 else if (pos[1] <= 0.15 && 1.20 <= pos[0])
526 template <
class Flu
idState,
class Context>
527 void initialFluidState_(FluidState& fs,
const Context& context,
528 unsigned spaceIdx,
unsigned timeIdx)
const
530 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
532 fs.setTemperature(293.0 );
536 if (isContaminated_(pos)) {
537 fs.setSaturation(waterPhaseIdx, 0.12);
538 fs.setSaturation(naplPhaseIdx, 0.07);
539 fs.setSaturation(gasPhaseIdx, 1 - 0.12 - 0.07);
543 Scalar pc[numPhases];
544 MaterialLaw::capillaryPressures(pc, matParams, fs);
545 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
546 fs.setPressure(phaseIdx, pw + (pc[phaseIdx] - pc[waterPhaseIdx]));
549 using MMPC = Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
550 typename FluidSystem::template ParameterCache<Scalar> paramCache;
551 MMPC::solve(fs, paramCache,
true,
true);
554 fs.setSaturation(waterPhaseIdx, 0.12);
555 fs.setSaturation(gasPhaseIdx, 1 - fs.saturation(waterPhaseIdx));
556 fs.setSaturation(naplPhaseIdx, 0);
560 Scalar pc[numPhases];
561 MaterialLaw::capillaryPressures(pc, matParams, fs);
562 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
563 fs.setPressure(phaseIdx, pw + (pc[phaseIdx] - pc[waterPhaseIdx]));
566 using MMPC = Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
567 typename FluidSystem::template ParameterCache<Scalar> paramCache;
568 MMPC::solve(fs, paramCache,
true,
true);
571 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
572 fs.setMoleFraction(phaseIdx, NAPLIdx, 0.0);
574 if (phaseIdx == naplPhaseIdx)
578 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
579 sumx += fs.moleFraction(phaseIdx, compIdx);
581 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
582 fs.setMoleFraction(phaseIdx, compIdx,
583 fs.moleFraction(phaseIdx, compIdx) / sumx);
588 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
590 Scalar lambdaGranite = 2.8;
593 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
594 fs.setTemperature(293.15);
595 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
596 fs.setPressure(phaseIdx, 1.0135e5);
599 typename FluidSystem::template ParameterCache<Scalar> paramCache;
600 paramCache.updateAll(fs);
601 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
602 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
603 fs.setDensity(phaseIdx, rho);
606 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
607 Scalar lambdaSaturated;
608 if (FluidSystem::isLiquid(phaseIdx)) {
609 Scalar lambdaFluid = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
611 std::pow(lambdaGranite, (1 - poro))
613 std::pow(lambdaFluid, poro);
616 lambdaSaturated = std::pow(lambdaGranite, (1 - poro));
618 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
619 if (!FluidSystem::isLiquid(phaseIdx))
620 params.setVacuumLambda(lambdaSaturated);
624 void initInjectFluidState_()
626 injectFluidState_.setTemperature(383.0);
627 injectFluidState_.setPressure(gasPhaseIdx, 1e5);
628 injectFluidState_.setSaturation(gasPhaseIdx, 1.0);
630 Scalar xgH2O = 0.417;
631 injectFluidState_.setMoleFraction(gasPhaseIdx, H2OIdx, xgH2O);
632 injectFluidState_.setMoleFraction(gasPhaseIdx, airIdx, 1 - xgH2O);
633 injectFluidState_.setMoleFraction(gasPhaseIdx, NAPLIdx, 0.0);
636 typename FluidSystem::template ParameterCache<Scalar> paramCache;
637 paramCache.updatePhase(injectFluidState_, gasPhaseIdx);
639 Scalar h = FluidSystem::enthalpy(injectFluidState_, paramCache, gasPhaseIdx);
640 injectFluidState_.setEnthalpy(gasPhaseIdx, h);
646 Scalar finePorosity_;
647 Scalar coarsePorosity_;
649 MaterialLawParams fineMaterialParams_;
650 MaterialLawParams coarseMaterialParams_;
652 ThermalConductionLawParams thermalCondParams_;
653 SolidEnergyLawParams solidEnergyLawParams_;
655 Opm::CompositionalFluidState<Scalar, FluidSystem> injectFluidState_;
Non-isothermal three-phase gas injection problem where a hot gas is injected into a unsaturated porou...
Definition: cuvetteproblem.hh:180
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:422
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:365
CuvetteProblem(Simulator &simulator)
Definition: cuvetteproblem.hh:220
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:471
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:378
void finishInit()
Definition: cuvetteproblem.hh:228
const ThermalConductionLawParams & thermalConductionParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: cuvetteproblem.hh:406
std::string name() const
Definition: cuvetteproblem.hh:323
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:391
void endTimeStep()
Definition: cuvetteproblem.hh:329
bool shouldWriteRestartFile() const
Definition: cuvetteproblem.hh:317
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: cuvetteproblem.hh:356
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: cuvetteproblem.hh:489
Definition: cuvetteproblem.hh:64