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>;
94 using type = Opm::ThreePhaseParkerVanGenuchten<Traits>;
98template<
class TypeTag>
100{
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
103template<
class TypeTag>
112 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
145template <
class TypeTag>
166 enum { numPhases = FluidSystem::numPhases };
167 enum { numComponents = FluidSystem::numComponents };
168 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
169 enum { naplPhaseIdx = FluidSystem::naplPhaseIdx };
170 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
171 enum { H2OIdx = FluidSystem::H2OIdx };
172 enum { airIdx = FluidSystem::airIdx };
173 enum { NAPLIdx = FluidSystem::NAPLIdx };
174 enum { conti0EqIdx = Indices::conti0EqIdx };
177 enum { dimWorld = GridView::dimensionworld };
179 using CoordScalar =
typename GridView::ctype;
180 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
181 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
188 : ParentType(simulator)
197 ParentType::finishInit();
199 if (Opm::Valgrind::IsRunning())
200 FluidSystem::init(283.15, 500.0, 20,
203 FluidSystem::init(283.15, 500.0, 200,
207 fineK_ = this->toDimMatrix_(6.28e-12);
208 coarseK_ = this->toDimMatrix_(9.14e-10);
211 finePorosity_ = 0.42;
212 coarsePorosity_ = 0.42;
217 fineMaterialParams_.setVgAlpha(0.0005);
218 coarseMaterialParams_.setVgAlpha(0.005);
219 fineMaterialParams_.setVgN(4.0);
220 coarseMaterialParams_.setVgN(4.0);
222 coarseMaterialParams_.setkrRegardsSnr(
true);
223 fineMaterialParams_.setkrRegardsSnr(
true);
226 fineMaterialParams_.setSwr(0.1201);
227 fineMaterialParams_.setSwrx(0.1201);
228 fineMaterialParams_.setSnr(0.0701);
229 fineMaterialParams_.setSgr(0.0101);
230 coarseMaterialParams_.setSwr(0.1201);
231 coarseMaterialParams_.setSwrx(0.1201);
232 coarseMaterialParams_.setSnr(0.0701);
233 coarseMaterialParams_.setSgr(0.0101);
236 fineMaterialParams_.setPcMinSat(gasPhaseIdx, 0);
237 fineMaterialParams_.setPcMaxSat(gasPhaseIdx, 0);
238 fineMaterialParams_.setPcMinSat(naplPhaseIdx, 0);
239 fineMaterialParams_.setPcMaxSat(naplPhaseIdx, -1000);
240 fineMaterialParams_.setPcMinSat(waterPhaseIdx, 0);
241 fineMaterialParams_.setPcMaxSat(waterPhaseIdx, -10000);
243 coarseMaterialParams_.setPcMinSat(gasPhaseIdx, 0);
244 coarseMaterialParams_.setPcMaxSat(gasPhaseIdx, 0);
245 coarseMaterialParams_.setPcMinSat(naplPhaseIdx, 0);
246 coarseMaterialParams_.setPcMaxSat(naplPhaseIdx, -100);
247 coarseMaterialParams_.setPcMinSat(waterPhaseIdx, 0);
248 coarseMaterialParams_.setPcMaxSat(waterPhaseIdx, -1000);
251 fineMaterialParams_.setResidSat(waterPhaseIdx, 0.1201);
252 fineMaterialParams_.setResidSat(naplPhaseIdx, 0.0701);
253 fineMaterialParams_.setResidSat(gasPhaseIdx, 0.0101);
255 coarseMaterialParams_.setResidSat(waterPhaseIdx, 0.1201);
256 coarseMaterialParams_.setResidSat(naplPhaseIdx, 0.0701);
257 coarseMaterialParams_.setResidSat(gasPhaseIdx, 0.0101);
260 fineMaterialParams_.finalize();
261 coarseMaterialParams_.finalize();
264 computeThermalCondParams_(thermalCondParams_, finePorosity_);
267 solidEnergyLawParams_.setSolidHeatCapacity(790.0
269 solidEnergyLawParams_.finalize();
271 initInjectFluidState_();
279 ParentType::registerParameters();
281 Parameters::SetDefault<Parameters::GridFile>(
"./data/cuvette_11x4.dgf");
282 Parameters::SetDefault<Parameters::EndTime<Scalar>>(100.0);
283 Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(1.0);
284 Parameters::SetDefault<Parameters::MaxTimeStepSize<Scalar>>(600.0);
285 Parameters::SetDefault<Parameters::EnableGravity>(
true);
305 {
return std::string(
"cuvette_") + Model::name(); }
313 this->model().checkConservativeness();
317 this->model().globalStorage(storage);
320 if (this->gridView().comm().rank() == 0) {
321 std::cout <<
"Storage: " << storage << std::endl << std::flush;
336 template <
class Context>
345 template <
class Context>
347 unsigned timeIdx)
const
349 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
350 if (isFineMaterial_(pos))
358 template <
class Context>
359 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
361 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
362 if (isFineMaterial_(pos))
363 return finePorosity_;
365 return coarsePorosity_;
371 template <
class Context>
373 unsigned spaceIdx,
unsigned timeIdx)
const
375 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
376 if (isFineMaterial_(pos))
377 return fineMaterialParams_;
379 return coarseMaterialParams_;
385 template <
class Context>
386 const ThermalConductionLawParams &
390 {
return thermalCondParams_; }
402 template <
class Context>
403 void boundary(BoundaryRateVector& values,
const Context& context,
404 unsigned spaceIdx,
unsigned timeIdx)
const
406 const auto& pos = context.pos(spaceIdx, timeIdx);
408 if (onRightBoundary_(pos)) {
409 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
411 initialFluidState_(fs, context, spaceIdx, timeIdx);
413 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
416 else if (onLeftBoundary_(pos)) {
418 RateVector molarRate;
422 Scalar molarInjectionRate = 0.3435;
423 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
424 molarRate[conti0EqIdx + compIdx] =
426 * injectFluidState_.moleFraction(gasPhaseIdx, compIdx);
429 Scalar massInjectionRate =
431 * injectFluidState_.averageMolarMass(gasPhaseIdx);
434 values.setMolarRate(molarRate);
435 values.setEnthalpyRate(-injectFluidState_.enthalpy(gasPhaseIdx) * massInjectionRate);
451 template <
class Context>
452 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
453 unsigned timeIdx)
const
455 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
457 initialFluidState_(fs, context, spaceIdx, timeIdx);
460 values.assignMassConservative(fs, matParams,
false);
469 template <
class Context>
474 { rate = Scalar(0.0); }
479 bool onLeftBoundary_(
const GlobalPosition& pos)
const
480 {
return pos[0] < eps_; }
482 bool onRightBoundary_(
const GlobalPosition& pos)
const
483 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
485 bool onLowerBoundary_(
const GlobalPosition& pos)
const
486 {
return pos[1] < eps_; }
488 bool onUpperBoundary_(
const GlobalPosition& pos)
const
489 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
491 bool isContaminated_(
const GlobalPosition& pos)
const
493 return (0.20 <= pos[0]) && (pos[0] <= 0.80) && (0.4 <= pos[1])
497 bool isFineMaterial_(
const GlobalPosition& pos)
const
499 if (0.13 <= pos[0] && 1.20 >= pos[0] && 0.32 <= pos[1] && pos[1] <= 0.57)
501 else if (pos[1] <= 0.15 && 1.20 <= pos[0])
507 template <
class Flu
idState,
class Context>
508 void initialFluidState_(FluidState& fs,
const Context& context,
509 unsigned spaceIdx,
unsigned timeIdx)
const
511 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
513 fs.setTemperature(293.0 );
517 if (isContaminated_(pos)) {
518 fs.setSaturation(waterPhaseIdx, 0.12);
519 fs.setSaturation(naplPhaseIdx, 0.07);
520 fs.setSaturation(gasPhaseIdx, 1 - 0.12 - 0.07);
524 Scalar pc[numPhases];
525 MaterialLaw::capillaryPressures(pc, matParams, fs);
526 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
527 fs.setPressure(phaseIdx, pw + (pc[phaseIdx] - pc[waterPhaseIdx]));
530 using MMPC = Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
531 typename FluidSystem::template ParameterCache<Scalar> paramCache;
532 MMPC::solve(fs, paramCache,
true,
true);
535 fs.setSaturation(waterPhaseIdx, 0.12);
536 fs.setSaturation(gasPhaseIdx, 1 - fs.saturation(waterPhaseIdx));
537 fs.setSaturation(naplPhaseIdx, 0);
541 Scalar pc[numPhases];
542 MaterialLaw::capillaryPressures(pc, matParams, fs);
543 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
544 fs.setPressure(phaseIdx, pw + (pc[phaseIdx] - pc[waterPhaseIdx]));
547 using MMPC = Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
548 typename FluidSystem::template ParameterCache<Scalar> paramCache;
549 MMPC::solve(fs, paramCache,
true,
true);
552 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
553 fs.setMoleFraction(phaseIdx, NAPLIdx, 0.0);
555 if (phaseIdx == naplPhaseIdx)
559 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
560 sumx += fs.moleFraction(phaseIdx, compIdx);
562 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
563 fs.setMoleFraction(phaseIdx, compIdx,
564 fs.moleFraction(phaseIdx, compIdx) / sumx);
569 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
571 Scalar lambdaGranite = 2.8;
574 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
575 fs.setTemperature(293.15);
576 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
577 fs.setPressure(phaseIdx, 1.0135e5);
580 typename FluidSystem::template ParameterCache<Scalar> paramCache;
581 paramCache.updateAll(fs);
582 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
583 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
584 fs.setDensity(phaseIdx, rho);
587 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
588 Scalar lambdaSaturated;
589 if (FluidSystem::isLiquid(phaseIdx)) {
590 Scalar lambdaFluid = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
592 std::pow(lambdaGranite, (1 - poro))
594 std::pow(lambdaFluid, poro);
597 lambdaSaturated = std::pow(lambdaGranite, (1 - poro));
599 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
600 if (!FluidSystem::isLiquid(phaseIdx))
601 params.setVacuumLambda(lambdaSaturated);
605 void initInjectFluidState_()
607 injectFluidState_.setTemperature(383.0);
608 injectFluidState_.setPressure(gasPhaseIdx, 1e5);
609 injectFluidState_.setSaturation(gasPhaseIdx, 1.0);
611 Scalar xgH2O = 0.417;
612 injectFluidState_.setMoleFraction(gasPhaseIdx, H2OIdx, xgH2O);
613 injectFluidState_.setMoleFraction(gasPhaseIdx, airIdx, 1 - xgH2O);
614 injectFluidState_.setMoleFraction(gasPhaseIdx, NAPLIdx, 0.0);
617 typename FluidSystem::template ParameterCache<Scalar> paramCache;
618 paramCache.updatePhase(injectFluidState_, gasPhaseIdx);
620 Scalar h = FluidSystem::enthalpy(injectFluidState_, paramCache, gasPhaseIdx);
621 injectFluidState_.setEnthalpy(gasPhaseIdx, h);
627 Scalar finePorosity_;
628 Scalar coarsePorosity_;
630 MaterialLawParams fineMaterialParams_;
631 MaterialLawParams coarseMaterialParams_;
633 ThermalConductionLawParams thermalCondParams_;
634 SolidEnergyLawParams solidEnergyLawParams_;
636 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:147
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: cuvetteproblem.hh:337
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: cuvetteproblem.hh:403
static void registerParameters()
Definition: cuvetteproblem.hh:277
CuvetteProblem(Simulator &simulator)
Definition: cuvetteproblem.hh:187
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:372
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: cuvetteproblem.hh:452
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:346
const ThermalConductionLawParams & thermalConductionParams(const Context &, unsigned, unsigned) const
Definition: cuvetteproblem.hh:387
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:359
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: cuvetteproblem.hh:195
std::string name() const
The problem name.
Definition: cuvetteproblem.hh:304
void endTimeStep()
Called by the simulator after each time integration.
Definition: cuvetteproblem.hh:310
bool shouldWriteRestartFile() const
Returns true if a restart file should be written to disk.
Definition: cuvetteproblem.hh:298
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:470
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
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:235
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:69
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:94
The material law which ought to be used (extracted from the spatial parameters)
Definition: multiphasebaseproperties.hh:51
The type of the problem.
Definition: fvbaseproperties.hh:81
Opm::ConstantSolidHeatCapLaw< GetPropType< TypeTag, Properties::Scalar > > type
Definition: cuvetteproblem.hh:100
The material law for the energy stored in the solid matrix.
Definition: multiphasebaseproperties.hh:57
Definition: cuvetteproblem.hh:63
Opm::SomertonThermalConductionLaw< FluidSystem, Scalar > type
Definition: cuvetteproblem.hh:112
The material law for thermal conduction.
Definition: multiphasebaseproperties.hh:63