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/heatconduction/Somerton.hpp> 39 #include <opm/material/constraintsolvers/MiscibleMultiPhaseComposition.hpp> 40 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp> 41 #include <opm/common/Valgrind.hpp> 42 #include <opm/common/Unused.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>
59 namespace Properties {
72 CuvetteBaseProblem, FluidSystem,
73 Opm::FluidSystems::H2OAirMesitylene<
typename GET_PROP_TYPE(TypeTag, Scalar)>);
82 SET_PROP(CuvetteBaseProblem, MaterialLaw)
86 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
88 typedef Opm::ThreePhaseMaterialTraits<
90 FluidSystem::waterPhaseIdx,
91 FluidSystem::naplPhaseIdx,
92 FluidSystem::gasPhaseIdx> Traits;
95 typedef Opm::ThreePhaseParkerVanGenuchten<Traits> type;
99 SET_PROP(CuvetteBaseProblem, HeatConductionLaw)
103 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
107 typedef Opm::Somerton<FluidSystem, Scalar> type;
117 SET_STRING_PROP(CuvetteBaseProblem, GridFile,
"./data/cuvette_11x4.dgf");
149 template <
class TypeTag>
150 class CuvetteProblem :
public GET_PROP_TYPE(TypeTag, BaseProblem)
152 typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
154 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
155 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
156 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
157 typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
158 typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLaw) HeatConductionLaw;
159 typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLawParams) HeatConductionLawParams;
160 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
161 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
162 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
163 typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
164 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
165 typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
166 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
169 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
170 enum { numPhases = FluidSystem::numPhases };
171 enum { numComponents = FluidSystem::numComponents };
172 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
173 enum { naplPhaseIdx = FluidSystem::naplPhaseIdx };
174 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
175 enum { H2OIdx = FluidSystem::H2OIdx };
176 enum { airIdx = FluidSystem::airIdx };
177 enum { NAPLIdx = FluidSystem::NAPLIdx };
178 enum { conti0EqIdx = Indices::conti0EqIdx };
181 enum { dimWorld = GridView::dimensionworld };
183 typedef typename GridView::ctype CoordScalar;
184 typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
185 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
192 : ParentType(simulator)
201 ParentType::finishInit();
203 if (Opm::Valgrind::IsRunning())
204 FluidSystem::init(283.15, 500.0, 20,
207 FluidSystem::init(283.15, 500.0, 200,
211 fineK_ = this->toDimMatrix_(6.28e-12);
212 coarseK_ = this->toDimMatrix_(9.14e-10);
215 finePorosity_ = 0.42;
216 coarsePorosity_ = 0.42;
221 fineMaterialParams_.setVgAlpha(0.0005);
222 coarseMaterialParams_.setVgAlpha(0.005);
223 fineMaterialParams_.setVgN(4.0);
224 coarseMaterialParams_.setVgN(4.0);
226 coarseMaterialParams_.setkrRegardsSnr(
true);
227 fineMaterialParams_.setkrRegardsSnr(
true);
230 fineMaterialParams_.setSwr(0.1201);
231 fineMaterialParams_.setSwrx(0.1201);
232 fineMaterialParams_.setSnr(0.0701);
233 fineMaterialParams_.setSgr(0.0101);
234 coarseMaterialParams_.setSwr(0.1201);
235 coarseMaterialParams_.setSwrx(0.1201);
236 coarseMaterialParams_.setSnr(0.0701);
237 coarseMaterialParams_.setSgr(0.0101);
240 fineMaterialParams_.setPcMinSat(gasPhaseIdx, 0);
241 fineMaterialParams_.setPcMaxSat(gasPhaseIdx, 0);
242 fineMaterialParams_.setPcMinSat(naplPhaseIdx, 0);
243 fineMaterialParams_.setPcMaxSat(naplPhaseIdx, -1000);
244 fineMaterialParams_.setPcMinSat(waterPhaseIdx, 0);
245 fineMaterialParams_.setPcMaxSat(waterPhaseIdx, -10000);
247 coarseMaterialParams_.setPcMinSat(gasPhaseIdx, 0);
248 coarseMaterialParams_.setPcMaxSat(gasPhaseIdx, 0);
249 coarseMaterialParams_.setPcMinSat(naplPhaseIdx, 0);
250 coarseMaterialParams_.setPcMaxSat(naplPhaseIdx, -100);
251 coarseMaterialParams_.setPcMinSat(waterPhaseIdx, 0);
252 coarseMaterialParams_.setPcMaxSat(waterPhaseIdx, -1000);
255 fineMaterialParams_.setResidSat(waterPhaseIdx, 0.1201);
256 fineMaterialParams_.setResidSat(naplPhaseIdx, 0.0701);
257 fineMaterialParams_.setResidSat(gasPhaseIdx, 0.0101);
259 coarseMaterialParams_.setResidSat(waterPhaseIdx, 0.1201);
260 coarseMaterialParams_.setResidSat(naplPhaseIdx, 0.0701);
261 coarseMaterialParams_.setResidSat(gasPhaseIdx, 0.0101);
264 fineMaterialParams_.finalize();
265 coarseMaterialParams_.finalize();
268 computeHeatCondParams_(heatCondParams_, finePorosity_);
270 initInjectFluidState_();
290 {
return std::string(
"cuvette_") + Model::name(); }
298 this->model().checkConservativeness();
302 this->model().globalStorage(storage);
305 if (this->gridView().comm().rank() == 0) {
306 std::cout <<
"Storage: " << storage << std::endl << std::flush;
321 template <
class Context>
323 unsigned spaceIdx OPM_UNUSED,
324 unsigned timeIdx OPM_UNUSED)
const 330 template <
class Context>
332 unsigned timeIdx)
const 334 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
335 if (isFineMaterial_(pos))
343 template <
class Context>
344 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const 346 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
347 if (isFineMaterial_(pos))
348 return finePorosity_;
350 return coarsePorosity_;
356 template <
class Context>
358 unsigned spaceIdx,
unsigned timeIdx)
const 360 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
361 if (isFineMaterial_(pos))
362 return fineMaterialParams_;
364 return coarseMaterialParams_;
370 template <
class Context>
371 const HeatConductionLawParams &
373 unsigned spaceIdx OPM_UNUSED,
374 unsigned timeIdx OPM_UNUSED)
const 375 {
return heatCondParams_; }
380 template <
class Context>
382 unsigned spaceIdx OPM_UNUSED,
383 unsigned timeIdx OPM_UNUSED)
const 399 template <
class Context>
400 void boundary(BoundaryRateVector& values,
const Context& context,
401 unsigned spaceIdx,
unsigned timeIdx)
const 403 const auto& pos = context.pos(spaceIdx, timeIdx);
405 if (onRightBoundary_(pos)) {
406 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
408 initialFluidState_(fs, context, spaceIdx, timeIdx);
410 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
413 else if (onLeftBoundary_(pos)) {
415 RateVector molarRate;
419 Scalar molarInjectionRate = 0.3435;
420 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
421 molarRate[conti0EqIdx + compIdx] =
423 * injectFluidState_.moleFraction(gasPhaseIdx, compIdx);
426 Scalar massInjectionRate =
428 * injectFluidState_.averageMolarMass(gasPhaseIdx);
431 values.setMolarRate(molarRate);
432 values.setEnthalpyRate(-injectFluidState_.enthalpy(gasPhaseIdx) * massInjectionRate);
448 template <
class Context>
449 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
450 unsigned timeIdx)
const 452 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
454 initialFluidState_(fs, context, spaceIdx, timeIdx);
457 values.assignMassConservative(fs, matParams,
false);
466 template <
class Context>
468 const Context& context OPM_UNUSED,
469 unsigned spaceIdx OPM_UNUSED,
470 unsigned timeIdx OPM_UNUSED)
const 471 { rate = Scalar(0.0); }
476 bool onLeftBoundary_(
const GlobalPosition& pos)
const 477 {
return pos[0] < eps_; }
479 bool onRightBoundary_(
const GlobalPosition& pos)
const 480 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
482 bool onLowerBoundary_(
const GlobalPosition& pos)
const 483 {
return pos[1] < eps_; }
485 bool onUpperBoundary_(
const GlobalPosition& pos)
const 486 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
488 bool isContaminated_(
const GlobalPosition& pos)
const 490 return (0.20 <= pos[0]) && (pos[0] <= 0.80) && (0.4 <= pos[1])
494 bool isFineMaterial_(
const GlobalPosition& pos)
const 496 if (0.13 <= pos[0] && 1.20 >= pos[0] && 0.32 <= pos[1] && pos[1] <= 0.57)
498 else if (pos[1] <= 0.15 && 1.20 <= pos[0])
504 template <
class Flu
idState,
class Context>
505 void initialFluidState_(FluidState& fs,
const Context& context,
506 unsigned spaceIdx,
unsigned timeIdx)
const 508 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
510 fs.setTemperature(293.0 );
514 if (isContaminated_(pos)) {
515 fs.setSaturation(waterPhaseIdx, 0.12);
516 fs.setSaturation(naplPhaseIdx, 0.07);
517 fs.setSaturation(gasPhaseIdx, 1 - 0.12 - 0.07);
521 Scalar pc[numPhases];
522 MaterialLaw::capillaryPressures(pc, matParams, fs);
523 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
524 fs.setPressure(phaseIdx, pw + (pc[phaseIdx] - pc[waterPhaseIdx]));
527 typedef Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MMPC;
528 typename FluidSystem::template ParameterCache<Scalar> paramCache;
529 MMPC::solve(fs, paramCache,
true,
true);
532 fs.setSaturation(waterPhaseIdx, 0.12);
533 fs.setSaturation(gasPhaseIdx, 1 - fs.saturation(waterPhaseIdx));
534 fs.setSaturation(naplPhaseIdx, 0);
538 Scalar pc[numPhases];
539 MaterialLaw::capillaryPressures(pc, matParams, fs);
540 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
541 fs.setPressure(phaseIdx, pw + (pc[phaseIdx] - pc[waterPhaseIdx]));
544 typedef Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MMPC;
545 typename FluidSystem::template ParameterCache<Scalar> paramCache;
546 MMPC::solve(fs, paramCache,
true,
true);
549 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
550 fs.setMoleFraction(phaseIdx, NAPLIdx, 0.0);
552 if (phaseIdx == naplPhaseIdx)
556 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
557 sumx += fs.moleFraction(phaseIdx, compIdx);
559 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
560 fs.setMoleFraction(phaseIdx, compIdx,
561 fs.moleFraction(phaseIdx, compIdx) / sumx);
566 void computeHeatCondParams_(HeatConductionLawParams& params, Scalar poro)
568 Scalar lambdaGranite = 2.8;
571 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
572 fs.setTemperature(293.15);
573 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
574 fs.setPressure(phaseIdx, 1.0135e5);
577 typename FluidSystem::template ParameterCache<Scalar> paramCache;
578 paramCache.updateAll(fs);
579 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
580 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
581 fs.setDensity(phaseIdx, rho);
584 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
585 Scalar lambdaSaturated;
586 if (FluidSystem::isLiquid(phaseIdx)) {
587 Scalar lambdaFluid = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
589 std::pow(lambdaGranite, (1 - poro))
591 std::pow(lambdaFluid, poro);
594 lambdaSaturated = std::pow(lambdaGranite, (1 - poro));
596 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
597 if (!FluidSystem::isLiquid(phaseIdx))
598 params.setVacuumLambda(lambdaSaturated);
602 void initInjectFluidState_()
604 injectFluidState_.setTemperature(383.0);
605 injectFluidState_.setPressure(gasPhaseIdx, 1e5);
606 injectFluidState_.setSaturation(gasPhaseIdx, 1.0);
608 Scalar xgH2O = 0.417;
609 injectFluidState_.setMoleFraction(gasPhaseIdx, H2OIdx, xgH2O);
610 injectFluidState_.setMoleFraction(gasPhaseIdx, airIdx, 1 - xgH2O);
611 injectFluidState_.setMoleFraction(gasPhaseIdx, NAPLIdx, 0.0);
614 typename FluidSystem::template ParameterCache<Scalar> paramCache;
615 paramCache.updatePhase(injectFluidState_, gasPhaseIdx);
617 Scalar h = FluidSystem::enthalpy(injectFluidState_, paramCache, gasPhaseIdx);
618 injectFluidState_.setEnthalpy(gasPhaseIdx, h);
624 Scalar finePorosity_;
625 Scalar coarsePorosity_;
627 MaterialLawParams fineMaterialParams_;
628 MaterialLawParams coarseMaterialParams_;
630 HeatConductionLawParams heatCondParams_;
632 Opm::CompositionalFluidState<Scalar, FluidSystem> injectFluidState_;
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: cuvetteproblem.hh:449
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
const HeatConductionLawParams & heatConductionParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: cuvetteproblem.hh:372
CuvetteProblem(Simulator &simulator)
Definition: cuvetteproblem.hh:191
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:344
Definition: baseauxiliarymodule.hh:37
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:331
void endTimeStep()
Called by the simulator after each time integration.
Definition: cuvetteproblem.hh:295
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
Non-isothermal three-phase gas injection problem where a hot gas is injected into a unsaturated porou...
Definition: cuvetteproblem.hh:55
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: cuvetteproblem.hh:322
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: cuvetteproblem.hh:199
bool shouldWriteRestartFile() const
Returns true if a restart file should be written to disk.
Definition: cuvetteproblem.hh:283
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: cuvetteproblem.hh:467
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: cuvetteproblem.hh:357
Scalar heatCapacitySolid(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: cuvetteproblem.hh:381
std::string name() const
The problem name.
Definition: cuvetteproblem.hh:289
Declares the properties required for the compositional multi-phase primary variable switching model...
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: cuvetteproblem.hh:400
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
#define SET_STRING_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant string value.
Definition: propertysystem.hh:416
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394