28#ifndef EWOMS_CUVETTE_PROBLEM_HH
29#define EWOMS_CUVETTE_PROBLEM_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>
54template <
class TypeTag>
67template<
class TypeTag>
68struct Grid<TypeTag, TTag::CuvetteBaseProblem> {
using type = Dune::YaspGrid<2>; };
71template<
class TypeTag>
75template<
class TypeTag>
77{
using type = Opm::H2OAirMesityleneFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
80template<
class TypeTag>
87 using Traits = Opm::ThreePhaseMaterialTraits<
89 FluidSystem::waterPhaseIdx,
90 FluidSystem::naplPhaseIdx,
91 FluidSystem::gasPhaseIdx,
96 using type = Opm::ThreePhaseParkerVanGenuchten<Traits>;
100template<
class TypeTag>
102{
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
105template<
class TypeTag>
114 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
147template <
class TypeTag>
168 enum { numPhases = FluidSystem::numPhases };
169 enum { numComponents = FluidSystem::numComponents };
170 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
171 enum { naplPhaseIdx = FluidSystem::naplPhaseIdx };
172 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
173 enum { H2OIdx = FluidSystem::H2OIdx };
174 enum { airIdx = FluidSystem::airIdx };
175 enum { NAPLIdx = FluidSystem::NAPLIdx };
176 enum { conti0EqIdx = Indices::conti0EqIdx };
179 enum { dimWorld = GridView::dimensionworld };
181 using CoordScalar =
typename GridView::ctype;
182 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
183 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
190 : ParentType(simulator)
199 ParentType::finishInit();
201 if (Opm::Valgrind::IsRunning())
202 FluidSystem::init(283.15, 500.0, 20,
205 FluidSystem::init(283.15, 500.0, 200,
209 fineK_ = this->toDimMatrix_(6.28e-12);
210 coarseK_ = this->toDimMatrix_(9.14e-10);
213 finePorosity_ = 0.42;
214 coarsePorosity_ = 0.42;
219 fineMaterialParams_.setVgAlpha(0.0005);
220 coarseMaterialParams_.setVgAlpha(0.005);
221 fineMaterialParams_.setVgN(4.0);
222 coarseMaterialParams_.setVgN(4.0);
224 coarseMaterialParams_.setkrRegardsSnr(
true);
225 fineMaterialParams_.setkrRegardsSnr(
true);
228 fineMaterialParams_.setSwr(0.1201);
229 fineMaterialParams_.setSwrx(0.1201);
230 fineMaterialParams_.setSnr(0.0701);
231 fineMaterialParams_.setSgr(0.0101);
232 coarseMaterialParams_.setSwr(0.1201);
233 coarseMaterialParams_.setSwrx(0.1201);
234 coarseMaterialParams_.setSnr(0.0701);
235 coarseMaterialParams_.setSgr(0.0101);
238 fineMaterialParams_.setPcMinSat(gasPhaseIdx, 0);
239 fineMaterialParams_.setPcMaxSat(gasPhaseIdx, 0);
240 fineMaterialParams_.setPcMinSat(naplPhaseIdx, 0);
241 fineMaterialParams_.setPcMaxSat(naplPhaseIdx, -1000);
242 fineMaterialParams_.setPcMinSat(waterPhaseIdx, 0);
243 fineMaterialParams_.setPcMaxSat(waterPhaseIdx, -10000);
245 coarseMaterialParams_.setPcMinSat(gasPhaseIdx, 0);
246 coarseMaterialParams_.setPcMaxSat(gasPhaseIdx, 0);
247 coarseMaterialParams_.setPcMinSat(naplPhaseIdx, 0);
248 coarseMaterialParams_.setPcMaxSat(naplPhaseIdx, -100);
249 coarseMaterialParams_.setPcMinSat(waterPhaseIdx, 0);
250 coarseMaterialParams_.setPcMaxSat(waterPhaseIdx, -1000);
253 fineMaterialParams_.setResidSat(waterPhaseIdx, 0.1201);
254 fineMaterialParams_.setResidSat(naplPhaseIdx, 0.0701);
255 fineMaterialParams_.setResidSat(gasPhaseIdx, 0.0101);
257 coarseMaterialParams_.setResidSat(waterPhaseIdx, 0.1201);
258 coarseMaterialParams_.setResidSat(naplPhaseIdx, 0.0701);
259 coarseMaterialParams_.setResidSat(gasPhaseIdx, 0.0101);
262 fineMaterialParams_.finalize();
263 coarseMaterialParams_.finalize();
266 computeThermalCondParams_(thermalCondParams_, finePorosity_);
269 solidEnergyLawParams_.setSolidHeatCapacity(790.0
271 solidEnergyLawParams_.finalize();
273 initInjectFluidState_();
281 ParentType::registerParameters();
283 Parameters::SetDefault<Parameters::GridFile>(
"./data/cuvette_11x4.dgf");
284 Parameters::SetDefault<Parameters::EndTime<Scalar>>(100.0);
285 Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(1.0);
286 Parameters::SetDefault<Parameters::MaxTimeStepSize<Scalar>>(600.0);
287 Parameters::SetDefault<Parameters::EnableGravity>(
true);
307 {
return std::string(
"cuvette_") + Model::name(); }
315 this->model().checkConservativeness();
319 this->model().globalStorage(storage);
322 if (this->gridView().comm().rank() == 0) {
323 std::cout <<
"Storage: " << storage << std::endl << std::flush;
338 template <
class Context>
347 template <
class Context>
349 unsigned timeIdx)
const
351 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
352 if (isFineMaterial_(pos))
360 template <
class Context>
361 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
363 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
364 if (isFineMaterial_(pos))
365 return finePorosity_;
367 return coarsePorosity_;
373 template <
class Context>
375 unsigned spaceIdx,
unsigned timeIdx)
const
377 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
378 if (isFineMaterial_(pos))
379 return fineMaterialParams_;
381 return coarseMaterialParams_;
387 template <
class Context>
388 const ThermalConductionLawParams &
392 {
return thermalCondParams_; }
404 template <
class Context>
405 void boundary(BoundaryRateVector& values,
const Context& context,
406 unsigned spaceIdx,
unsigned timeIdx)
const
408 const auto& pos = context.pos(spaceIdx, timeIdx);
410 if (onRightBoundary_(pos)) {
411 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
413 initialFluidState_(fs, context, spaceIdx, timeIdx);
415 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
418 else if (onLeftBoundary_(pos)) {
420 RateVector molarRate;
424 Scalar molarInjectionRate = 0.3435;
425 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
426 molarRate[conti0EqIdx + compIdx] =
428 * injectFluidState_.moleFraction(gasPhaseIdx, compIdx);
431 Scalar massInjectionRate =
433 * injectFluidState_.averageMolarMass(gasPhaseIdx);
436 values.setMolarRate(molarRate);
437 values.setEnthalpyRate(-injectFluidState_.enthalpy(gasPhaseIdx) * massInjectionRate);
453 template <
class Context>
454 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
455 unsigned timeIdx)
const
457 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
459 initialFluidState_(fs, context, spaceIdx, timeIdx);
462 values.assignMassConservative(fs, matParams,
false);
471 template <
class Context>
476 { rate = Scalar(0.0); }
481 bool onLeftBoundary_(
const GlobalPosition& pos)
const
482 {
return pos[0] < eps_; }
484 bool onRightBoundary_(
const GlobalPosition& pos)
const
485 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
487 bool onLowerBoundary_(
const GlobalPosition& pos)
const
488 {
return pos[1] < eps_; }
490 bool onUpperBoundary_(
const GlobalPosition& pos)
const
491 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
493 bool isContaminated_(
const GlobalPosition& pos)
const
495 return (0.20 <= pos[0]) && (pos[0] <= 0.80) && (0.4 <= pos[1])
499 bool isFineMaterial_(
const GlobalPosition& pos)
const
501 if (0.13 <= pos[0] && 1.20 >= pos[0] && 0.32 <= pos[1] && pos[1] <= 0.57)
503 else if (pos[1] <= 0.15 && 1.20 <= pos[0])
509 template <
class Flu
idState,
class Context>
510 void initialFluidState_(FluidState& fs,
const Context& context,
511 unsigned spaceIdx,
unsigned timeIdx)
const
513 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
515 fs.setTemperature(293.0 );
519 if (isContaminated_(pos)) {
520 fs.setSaturation(waterPhaseIdx, 0.12);
521 fs.setSaturation(naplPhaseIdx, 0.07);
522 fs.setSaturation(gasPhaseIdx, 1 - 0.12 - 0.07);
526 Scalar pc[numPhases];
527 MaterialLaw::capillaryPressures(pc, matParams, fs);
528 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
529 fs.setPressure(phaseIdx, pw + (pc[phaseIdx] - pc[waterPhaseIdx]));
532 using MMPC = Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
533 typename FluidSystem::template ParameterCache<Scalar> paramCache;
534 MMPC::solve(fs, paramCache,
true,
true);
537 fs.setSaturation(waterPhaseIdx, 0.12);
538 fs.setSaturation(gasPhaseIdx, 1 - fs.saturation(waterPhaseIdx));
539 fs.setSaturation(naplPhaseIdx, 0);
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 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
555 fs.setMoleFraction(phaseIdx, NAPLIdx, 0.0);
557 if (phaseIdx == naplPhaseIdx)
561 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
562 sumx += fs.moleFraction(phaseIdx, compIdx);
564 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
565 fs.setMoleFraction(phaseIdx, compIdx,
566 fs.moleFraction(phaseIdx, compIdx) / sumx);
571 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
573 Scalar lambdaGranite = 2.8;
576 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
577 fs.setTemperature(293.15);
578 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
579 fs.setPressure(phaseIdx, 1.0135e5);
582 typename FluidSystem::template ParameterCache<Scalar> paramCache;
583 paramCache.updateAll(fs);
584 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
585 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
586 fs.setDensity(phaseIdx, rho);
589 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
590 Scalar lambdaSaturated;
591 if (FluidSystem::isLiquid(phaseIdx)) {
592 Scalar lambdaFluid = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
594 std::pow(lambdaGranite, (1 - poro))
596 std::pow(lambdaFluid, poro);
599 lambdaSaturated = std::pow(lambdaGranite, (1 - poro));
601 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
602 if (!FluidSystem::isLiquid(phaseIdx))
603 params.setVacuumLambda(lambdaSaturated);
607 void initInjectFluidState_()
609 injectFluidState_.setTemperature(383.0);
610 injectFluidState_.setPressure(gasPhaseIdx, 1e5);
611 injectFluidState_.setSaturation(gasPhaseIdx, 1.0);
613 Scalar xgH2O = 0.417;
614 injectFluidState_.setMoleFraction(gasPhaseIdx, H2OIdx, xgH2O);
615 injectFluidState_.setMoleFraction(gasPhaseIdx, airIdx, 1 - xgH2O);
616 injectFluidState_.setMoleFraction(gasPhaseIdx, NAPLIdx, 0.0);
619 typename FluidSystem::template ParameterCache<Scalar> paramCache;
620 paramCache.updatePhase(injectFluidState_, gasPhaseIdx);
622 Scalar h = FluidSystem::enthalpy(injectFluidState_, paramCache, gasPhaseIdx);
623 injectFluidState_.setEnthalpy(gasPhaseIdx, h);
629 Scalar finePorosity_;
630 Scalar coarsePorosity_;
632 MaterialLawParams fineMaterialParams_;
633 MaterialLawParams coarseMaterialParams_;
635 ThermalConductionLawParams thermalCondParams_;
636 SolidEnergyLawParams solidEnergyLawParams_;
638 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:149
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: cuvetteproblem.hh:339
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: cuvetteproblem.hh:405
static void registerParameters()
Definition: cuvetteproblem.hh:279
CuvetteProblem(Simulator &simulator)
Definition: cuvetteproblem.hh:189
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:374
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: cuvetteproblem.hh:454
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:348
const ThermalConductionLawParams & thermalConductionParams(const Context &, unsigned, unsigned) const
Definition: cuvetteproblem.hh:389
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:361
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: cuvetteproblem.hh:197
std::string name() const
The problem name.
Definition: cuvetteproblem.hh:306
void endTimeStep()
Called by the simulator after each time integration.
Definition: cuvetteproblem.hh:312
bool shouldWriteRestartFile() const
Returns true if a restart file should be written to disk.
Definition: cuvetteproblem.hh:300
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: cuvetteproblem.hh:472
Definition: blackoilmodel.hh:79
Definition: blackoilbioeffectsmodules.hh:43
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:233
Declares the properties required for the compositional multi-phase primary variable switching model.
Opm::H2OAirMesityleneFluidSystem< GetPropType< TypeTag, Properties::Scalar > > type
Definition: cuvetteproblem.hh:77
The fluid systems including the information about the phases.
Definition: multiphasebaseproperties.hh:79
Dune::YaspGrid< 2 > type
Definition: cuvetteproblem.hh:68
The type of the DUNE grid.
Definition: basicproperties.hh:100
Opm::ThreePhaseParkerVanGenuchten< Traits > type
Definition: cuvetteproblem.hh:96
The material law which ought to be used (extracted from the spatial parameters)
Definition: multiphasebaseproperties.hh:55
The type of the problem.
Definition: fvbaseproperties.hh:81
Opm::ConstantSolidHeatCapLaw< GetPropType< TypeTag, Properties::Scalar > > type
Definition: cuvetteproblem.hh:102
The material law for the energy stored in the solid matrix.
Definition: multiphasebaseproperties.hh:63
Definition: cuvetteproblem.hh:63
Opm::SomertonThermalConductionLaw< FluidSystem, Scalar > type
Definition: cuvetteproblem.hh:114
The material law for thermal conduction.
Definition: multiphasebaseproperties.hh:71