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>
44 #include <dune/grid/yaspgrid.hh>
45 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
47 #include <dune/common/version.hh>
48 #include <dune/common/fvector.hh>
49 #include <dune/common/fmatrix.hh>
54 template <
class TypeTag>
58 namespace Opm::Properties {
67 template<
class TypeTag>
68 struct Grid<TypeTag, TTag::CuvetteBaseProblem> {
using type = Dune::YaspGrid<2>; };
71 template<
class TypeTag>
75 template<
class TypeTag>
76 struct FluidSystem<TypeTag, TTag::CuvetteBaseProblem>
77 {
using type = Opm::H2OAirMesityleneFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
80 template<
class TypeTag>
81 struct EnableGravity<TypeTag, TTag::CuvetteBaseProblem> {
static constexpr
bool value =
true; };
84 template<
class TypeTag>
85 struct MaxTimeStepSize<TypeTag, TTag::CuvetteBaseProblem>
87 using type = GetPropType<TypeTag, Scalar>;
88 static constexpr type value = 600.;
92 template<
class TypeTag>
93 struct MaterialLaw<TypeTag, TTag::CuvetteBaseProblem>
96 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
97 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
99 using Traits = Opm::ThreePhaseMaterialTraits<
101 FluidSystem::waterPhaseIdx,
102 FluidSystem::naplPhaseIdx,
103 FluidSystem::gasPhaseIdx>;
106 using type = Opm::ThreePhaseParkerVanGenuchten<Traits>;
110 template<
class TypeTag>
111 struct SolidEnergyLaw<TypeTag, TTag::CuvetteBaseProblem>
112 {
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
115 template<
class TypeTag>
116 struct ThermalConductionLaw<TypeTag, TTag::CuvetteBaseProblem>
119 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
120 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
124 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
128 template<
class TypeTag>
129 struct EndTime<TypeTag, TTag::CuvetteBaseProblem>
131 using type = GetPropType<TypeTag, Scalar>;
132 static constexpr type value = 180;
136 template<
class TypeTag>
137 struct InitialTimeStepSize<TypeTag, TTag::CuvetteBaseProblem>
139 using type = GetPropType<TypeTag, Scalar>;
140 static constexpr type value = 1;
144 template<
class TypeTag>
145 struct GridFile<TypeTag, TTag::CuvetteBaseProblem> {
static constexpr
auto value =
"./data/cuvette_11x4.dgf"; };
177 template <
class TypeTag>
180 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
182 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
183 using GridView = GetPropType<TypeTag, Properties::GridView>;
184 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
185 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
186 using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
187 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
188 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
189 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
190 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
191 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
192 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
193 using Model = GetPropType<TypeTag, Properties::Model>;
194 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
197 using Indices = GetPropType<TypeTag, Properties::Indices>;
198 enum { numPhases = FluidSystem::numPhases };
199 enum { numComponents = FluidSystem::numComponents };
200 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
201 enum { naplPhaseIdx = FluidSystem::naplPhaseIdx };
202 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
203 enum { H2OIdx = FluidSystem::H2OIdx };
204 enum { airIdx = FluidSystem::airIdx };
205 enum { NAPLIdx = FluidSystem::NAPLIdx };
206 enum { conti0EqIdx = Indices::conti0EqIdx };
209 enum { dimWorld = GridView::dimensionworld };
211 using CoordScalar =
typename GridView::ctype;
212 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
213 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
220 : ParentType(simulator)
229 ParentType::finishInit();
231 if (Opm::Valgrind::IsRunning())
232 FluidSystem::init(283.15, 500.0, 20,
235 FluidSystem::init(283.15, 500.0, 200,
239 fineK_ = this->toDimMatrix_(6.28e-12);
240 coarseK_ = this->toDimMatrix_(9.14e-10);
243 finePorosity_ = 0.42;
244 coarsePorosity_ = 0.42;
249 fineMaterialParams_.setVgAlpha(0.0005);
250 coarseMaterialParams_.setVgAlpha(0.005);
251 fineMaterialParams_.setVgN(4.0);
252 coarseMaterialParams_.setVgN(4.0);
254 coarseMaterialParams_.setkrRegardsSnr(
true);
255 fineMaterialParams_.setkrRegardsSnr(
true);
258 fineMaterialParams_.setSwr(0.1201);
259 fineMaterialParams_.setSwrx(0.1201);
260 fineMaterialParams_.setSnr(0.0701);
261 fineMaterialParams_.setSgr(0.0101);
262 coarseMaterialParams_.setSwr(0.1201);
263 coarseMaterialParams_.setSwrx(0.1201);
264 coarseMaterialParams_.setSnr(0.0701);
265 coarseMaterialParams_.setSgr(0.0101);
268 fineMaterialParams_.setPcMinSat(gasPhaseIdx, 0);
269 fineMaterialParams_.setPcMaxSat(gasPhaseIdx, 0);
270 fineMaterialParams_.setPcMinSat(naplPhaseIdx, 0);
271 fineMaterialParams_.setPcMaxSat(naplPhaseIdx, -1000);
272 fineMaterialParams_.setPcMinSat(waterPhaseIdx, 0);
273 fineMaterialParams_.setPcMaxSat(waterPhaseIdx, -10000);
275 coarseMaterialParams_.setPcMinSat(gasPhaseIdx, 0);
276 coarseMaterialParams_.setPcMaxSat(gasPhaseIdx, 0);
277 coarseMaterialParams_.setPcMinSat(naplPhaseIdx, 0);
278 coarseMaterialParams_.setPcMaxSat(naplPhaseIdx, -100);
279 coarseMaterialParams_.setPcMinSat(waterPhaseIdx, 0);
280 coarseMaterialParams_.setPcMaxSat(waterPhaseIdx, -1000);
283 fineMaterialParams_.setResidSat(waterPhaseIdx, 0.1201);
284 fineMaterialParams_.setResidSat(naplPhaseIdx, 0.0701);
285 fineMaterialParams_.setResidSat(gasPhaseIdx, 0.0101);
287 coarseMaterialParams_.setResidSat(waterPhaseIdx, 0.1201);
288 coarseMaterialParams_.setResidSat(naplPhaseIdx, 0.0701);
289 coarseMaterialParams_.setResidSat(gasPhaseIdx, 0.0101);
292 fineMaterialParams_.finalize();
293 coarseMaterialParams_.finalize();
296 computeThermalCondParams_(thermalCondParams_, finePorosity_);
299 solidEnergyLawParams_.setSolidHeatCapacity(790.0
301 solidEnergyLawParams_.finalize();
303 initInjectFluidState_();
323 {
return std::string(
"cuvette_") + Model::name(); }
331 this->model().checkConservativeness();
335 this->model().globalStorage(storage);
338 if (this->gridView().comm().rank() == 0) {
339 std::cout <<
"Storage: " << storage << std::endl << std::flush;
354 template <
class Context>
363 template <
class Context>
365 unsigned timeIdx)
const
367 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
368 if (isFineMaterial_(pos))
376 template <
class Context>
377 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
379 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
380 if (isFineMaterial_(pos))
381 return finePorosity_;
383 return coarsePorosity_;
389 template <
class Context>
391 unsigned spaceIdx,
unsigned timeIdx)
const
393 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
394 if (isFineMaterial_(pos))
395 return fineMaterialParams_;
397 return coarseMaterialParams_;
403 template <
class Context>
404 const ThermalConductionLawParams &
408 {
return thermalCondParams_; }
420 template <
class Context>
421 void boundary(BoundaryRateVector& values,
const Context& context,
422 unsigned spaceIdx,
unsigned timeIdx)
const
424 const auto& pos = context.pos(spaceIdx, timeIdx);
426 if (onRightBoundary_(pos)) {
427 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
429 initialFluidState_(fs, context, spaceIdx, timeIdx);
431 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
434 else if (onLeftBoundary_(pos)) {
436 RateVector molarRate;
440 Scalar molarInjectionRate = 0.3435;
441 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
442 molarRate[conti0EqIdx + compIdx] =
444 * injectFluidState_.moleFraction(gasPhaseIdx, compIdx);
447 Scalar massInjectionRate =
449 * injectFluidState_.averageMolarMass(gasPhaseIdx);
452 values.setMolarRate(molarRate);
453 values.setEnthalpyRate(-injectFluidState_.enthalpy(gasPhaseIdx) * massInjectionRate);
469 template <
class Context>
470 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
471 unsigned timeIdx)
const
473 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
475 initialFluidState_(fs, context, spaceIdx, timeIdx);
478 values.assignMassConservative(fs, matParams,
false);
487 template <
class Context>
492 { rate = Scalar(0.0); }
497 bool onLeftBoundary_(
const GlobalPosition& pos)
const
498 {
return pos[0] < eps_; }
500 bool onRightBoundary_(
const GlobalPosition& pos)
const
501 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
503 bool onLowerBoundary_(
const GlobalPosition& pos)
const
504 {
return pos[1] < eps_; }
506 bool onUpperBoundary_(
const GlobalPosition& pos)
const
507 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
509 bool isContaminated_(
const GlobalPosition& pos)
const
511 return (0.20 <= pos[0]) && (pos[0] <= 0.80) && (0.4 <= pos[1])
515 bool isFineMaterial_(
const GlobalPosition& pos)
const
517 if (0.13 <= pos[0] && 1.20 >= pos[0] && 0.32 <= pos[1] && pos[1] <= 0.57)
519 else if (pos[1] <= 0.15 && 1.20 <= pos[0])
525 template <
class Flu
idState,
class Context>
526 void initialFluidState_(FluidState& fs,
const Context& context,
527 unsigned spaceIdx,
unsigned timeIdx)
const
529 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
531 fs.setTemperature(293.0 );
535 if (isContaminated_(pos)) {
536 fs.setSaturation(waterPhaseIdx, 0.12);
537 fs.setSaturation(naplPhaseIdx, 0.07);
538 fs.setSaturation(gasPhaseIdx, 1 - 0.12 - 0.07);
542 Scalar pc[numPhases];
543 MaterialLaw::capillaryPressures(pc, matParams, fs);
544 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
545 fs.setPressure(phaseIdx, pw + (pc[phaseIdx] - pc[waterPhaseIdx]));
548 using MMPC = Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
549 typename FluidSystem::template ParameterCache<Scalar> paramCache;
550 MMPC::solve(fs, paramCache,
true,
true);
553 fs.setSaturation(waterPhaseIdx, 0.12);
554 fs.setSaturation(gasPhaseIdx, 1 - fs.saturation(waterPhaseIdx));
555 fs.setSaturation(naplPhaseIdx, 0);
559 Scalar pc[numPhases];
560 MaterialLaw::capillaryPressures(pc, matParams, fs);
561 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
562 fs.setPressure(phaseIdx, pw + (pc[phaseIdx] - pc[waterPhaseIdx]));
565 using MMPC = Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
566 typename FluidSystem::template ParameterCache<Scalar> paramCache;
567 MMPC::solve(fs, paramCache,
true,
true);
570 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
571 fs.setMoleFraction(phaseIdx, NAPLIdx, 0.0);
573 if (phaseIdx == naplPhaseIdx)
577 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
578 sumx += fs.moleFraction(phaseIdx, compIdx);
580 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
581 fs.setMoleFraction(phaseIdx, compIdx,
582 fs.moleFraction(phaseIdx, compIdx) / sumx);
587 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
589 Scalar lambdaGranite = 2.8;
592 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
593 fs.setTemperature(293.15);
594 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
595 fs.setPressure(phaseIdx, 1.0135e5);
598 typename FluidSystem::template ParameterCache<Scalar> paramCache;
599 paramCache.updateAll(fs);
600 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
601 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
602 fs.setDensity(phaseIdx, rho);
605 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
606 Scalar lambdaSaturated;
607 if (FluidSystem::isLiquid(phaseIdx)) {
608 Scalar lambdaFluid = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
610 std::pow(lambdaGranite, (1 - poro))
612 std::pow(lambdaFluid, poro);
615 lambdaSaturated = std::pow(lambdaGranite, (1 - poro));
617 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
618 if (!FluidSystem::isLiquid(phaseIdx))
619 params.setVacuumLambda(lambdaSaturated);
623 void initInjectFluidState_()
625 injectFluidState_.setTemperature(383.0);
626 injectFluidState_.setPressure(gasPhaseIdx, 1e5);
627 injectFluidState_.setSaturation(gasPhaseIdx, 1.0);
629 Scalar xgH2O = 0.417;
630 injectFluidState_.setMoleFraction(gasPhaseIdx, H2OIdx, xgH2O);
631 injectFluidState_.setMoleFraction(gasPhaseIdx, airIdx, 1 - xgH2O);
632 injectFluidState_.setMoleFraction(gasPhaseIdx, NAPLIdx, 0.0);
635 typename FluidSystem::template ParameterCache<Scalar> paramCache;
636 paramCache.updatePhase(injectFluidState_, gasPhaseIdx);
638 Scalar h = FluidSystem::enthalpy(injectFluidState_, paramCache, gasPhaseIdx);
639 injectFluidState_.setEnthalpy(gasPhaseIdx, h);
645 Scalar finePorosity_;
646 Scalar coarsePorosity_;
648 MaterialLawParams fineMaterialParams_;
649 MaterialLawParams coarseMaterialParams_;
651 ThermalConductionLawParams thermalCondParams_;
652 SolidEnergyLawParams solidEnergyLawParams_;
654 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:179
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: cuvetteproblem.hh:355
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:421
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:364
CuvetteProblem(Simulator &simulator)
Definition: cuvetteproblem.hh:219
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:470
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:377
void finishInit()
Definition: cuvetteproblem.hh:227
std::string name() const
Definition: cuvetteproblem.hh:322
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:390
void endTimeStep()
Definition: cuvetteproblem.hh:328
bool shouldWriteRestartFile() const
Definition: cuvetteproblem.hh:316
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: cuvetteproblem.hh:488
const ThermalConductionLawParams & thermalConductionParams(const Context &, unsigned, unsigned) const
Definition: cuvetteproblem.hh:405
Definition: cuvetteproblem.hh:63