261 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
263 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
264 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
265 using GridView = GetPropType<TypeTag, Properties::GridView>;
266 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
268 enum { dim = GridView::dimension };
269 enum { dimWorld = GridView::dimensionworld };
272 using Indices = GetPropType<TypeTag, Properties::Indices>;
273 enum { numPhases = FluidSystem::numPhases };
274 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
275 enum { liquidPhaseIdx = FluidSystem::liquidPhaseIdx };
276 enum { CO2Idx = FluidSystem::CO2Idx };
277 enum { BrineIdx = FluidSystem::BrineIdx };
278 enum { conti0EqIdx = Indices::conti0EqIdx };
279 enum { contiCO2EqIdx = conti0EqIdx + CO2Idx };
281 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
282 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
283 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
284 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
285 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
286 using Model = GetPropType<TypeTag, Properties::Model>;
287 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
288 using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
289 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
290 using ThermalConductionLawParams =
typename ThermalConductionLaw::Params;
292 using Toolbox = Opm::MathToolbox<Evaluation>;
293 using CoordScalar =
typename GridView::ctype;
294 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
295 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
302 : ParentType(simulator)
310 ParentType::finishInit();
314 temperatureLow_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureLow);
315 temperatureHigh_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureHigh);
316 nTemperature_ = EWOMS_GET_PARAM(TypeTag,
unsigned, FluidSystemNumTemperature);
318 pressureLow_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemPressureLow);
319 pressureHigh_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemPressureHigh);
320 nPressure_ = EWOMS_GET_PARAM(TypeTag,
unsigned, FluidSystemNumPressure);
322 maxDepth_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxDepth);
323 temperature_ = EWOMS_GET_PARAM(TypeTag, Scalar, Temperature);
327 FluidSystem::init(temperatureLow_,
334 fineLayerBottom_ = 22.0;
337 fineK_ = this->toDimMatrix_(1e-13);
338 coarseK_ = this->toDimMatrix_(1e-12);
342 coarsePorosity_ = 0.3;
345 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
346 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
347 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
348 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
351 fineMaterialParams_.setEntryPressure(1e4);
352 coarseMaterialParams_.setEntryPressure(5e3);
353 fineMaterialParams_.setLambda(2.0);
354 coarseMaterialParams_.setLambda(2.0);
356 fineMaterialParams_.finalize();
357 coarseMaterialParams_.finalize();
360 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
361 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
364 solidEnergyLawParams_.setSolidHeatCapacity(790.0
366 solidEnergyLawParams_.finalize();
374 ParentType::registerParameters();
376 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemTemperatureLow,
377 "The lower temperature [K] for tabulation of the "
379 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemTemperatureHigh,
380 "The upper temperature [K] for tabulation of the "
382 EWOMS_REGISTER_PARAM(TypeTag,
unsigned, FluidSystemNumTemperature,
383 "The number of intervals between the lower and "
384 "upper temperature");
386 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemPressureLow,
387 "The lower pressure [Pa] for tabulation of the "
389 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemPressureHigh,
390 "The upper pressure [Pa] for tabulation of the "
392 EWOMS_REGISTER_PARAM(TypeTag,
unsigned, FluidSystemNumPressure,
393 "The number of intervals between the lower and "
396 EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature,
397 "The temperature [K] in the reservoir");
398 EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxDepth,
399 "The maximum depth [m] of the reservoir");
400 EWOMS_REGISTER_PARAM(TypeTag, std::string, SimulationName,
401 "The name of the simulation used for the output "
415 std::ostringstream oss;
416 oss << EWOMS_GET_PARAM(TypeTag, std::string, SimulationName)
417 <<
"_" << Model::name();
418 if (getPropValue<TypeTag, Properties::EnableEnergy>())
420 oss <<
"_" << Model::discretizationName();
430 Scalar tol = this->model().newtonMethod().tolerance()*1e5;
431 this->model().checkConservativeness(tol);
434 PrimaryVariables storageL, storageG;
435 this->model().globalPhaseStorage(storageL, 0);
436 this->model().globalPhaseStorage(storageG, 1);
439 if (this->gridView().comm().rank() == 0) {
440 std::cout <<
"Storage: liquid=[" << storageL <<
"]"
441 <<
" gas=[" << storageG <<
"]\n" << std::flush;
449 template <
class Context>
450 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
452 const auto& pos = context.pos(spaceIdx, timeIdx);
453 if (inHighTemperatureRegion_(pos))
454 return temperature_ + 100;
461 template <
class Context>
463 unsigned timeIdx)
const
465 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
466 if (isFineMaterial_(pos))
474 template <
class Context>
475 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
477 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
478 if (isFineMaterial_(pos))
479 return finePorosity_;
480 return coarsePorosity_;
486 template <
class Context>
488 unsigned spaceIdx,
unsigned timeIdx)
const
490 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
491 if (isFineMaterial_(pos))
492 return fineMaterialParams_;
493 return coarseMaterialParams_;
501 template <
class Context>
502 const SolidEnergyLawParams&
506 {
return solidEnergyLawParams_; }
511 template <
class Context>
512 const ThermalConductionLawParams &
515 unsigned timeIdx)
const
517 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
518 if (isFineMaterial_(pos))
519 return fineThermalCondParams_;
520 return coarseThermalCondParams_;
533 template <
class Context>
534 void boundary(BoundaryRateVector& values,
const Context& context,
535 unsigned spaceIdx,
unsigned timeIdx)
const
537 const auto& pos = context.pos(spaceIdx, timeIdx);
538 if (onLeftBoundary_(pos)) {
539 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
540 initialFluidState_(fs, context, spaceIdx, timeIdx);
544 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
546 else if (onInlet_(pos)) {
547 RateVector massRate(0.0);
548 massRate[contiCO2EqIdx] = -1e-3;
550 using FluidState = Opm::ImmiscibleFluidState<Scalar, FluidSystem>;
552 fs.setSaturation(gasPhaseIdx, 1.0);
554 context.intensiveQuantities(spaceIdx, timeIdx).fluidState().pressure(gasPhaseIdx);
555 fs.setPressure(gasPhaseIdx, Toolbox::value(pg));
556 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
558 typename FluidSystem::template ParameterCache<Scalar> paramCache;
559 paramCache.updatePhase(fs, gasPhaseIdx);
560 Scalar h = FluidSystem::template enthalpy<FluidState, Scalar>(fs, paramCache, gasPhaseIdx);
563 values.setMassRate(massRate);
564 values.setEnthalpyRate(massRate[contiCO2EqIdx] * h);
581 template <
class Context>
582 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
583 unsigned timeIdx)
const
585 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
586 initialFluidState_(fs, context, spaceIdx, timeIdx);
591 values.assignNaive(fs);
600 template <
class Context>
605 { rate = Scalar(0.0); }
610 template <
class Context,
class Flu
idState>
611 void initialFluidState_(FluidState& fs,
612 const Context& context,
614 unsigned timeIdx)
const
616 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
621 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
626 fs.setSaturation(FluidSystem::liquidPhaseIdx, 1.0);
627 fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
632 Scalar densityL = FluidSystem::Brine::liquidDensity(temperature_, Scalar(1e5));
633 Scalar depth = maxDepth_ - pos[dim - 1];
634 Scalar pl = 1e5 - densityL * this->gravity()[dim - 1] * depth;
636 Scalar pC[numPhases];
638 MaterialLaw::capillaryPressures(pC, matParams, fs);
640 fs.setPressure(liquidPhaseIdx, pl + (pC[liquidPhaseIdx] - pC[liquidPhaseIdx]));
641 fs.setPressure(gasPhaseIdx, pl + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
646 fs.setMoleFraction(liquidPhaseIdx, CO2Idx, 0.005);
647 fs.setMoleFraction(liquidPhaseIdx, BrineIdx,
648 1.0 - fs.moleFraction(liquidPhaseIdx, CO2Idx));
650 typename FluidSystem::template ParameterCache<Scalar> paramCache;
651 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
652 CFRP::solve(fs, paramCache,
658 bool onLeftBoundary_(
const GlobalPosition& pos)
const
659 {
return pos[0] < eps_; }
661 bool onRightBoundary_(
const GlobalPosition& pos)
const
662 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
664 bool onInlet_(
const GlobalPosition& pos)
const
665 {
return onRightBoundary_(pos) && (5 < pos[1]) && (pos[1] < 15); }
667 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const
668 {
return (pos[0] > 20) && (pos[0] < 30) && (pos[1] > 5) && (pos[1] < 35); }
670 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
672 Scalar lambdaWater = 0.6;
673 Scalar lambdaGranite = 2.8;
675 Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
676 * std::pow(lambdaWater, poro);
677 Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
679 params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
680 params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
681 params.setVacuumLambda(lambdaDry);
684 bool isFineMaterial_(
const GlobalPosition& pos)
const
685 {
return pos[dim - 1] > fineLayerBottom_; }
689 Scalar fineLayerBottom_;
691 Scalar finePorosity_;
692 Scalar coarsePorosity_;
694 MaterialLawParams fineMaterialParams_;
695 MaterialLawParams coarseMaterialParams_;
697 ThermalConductionLawParams fineThermalCondParams_;
698 ThermalConductionLawParams coarseThermalCondParams_;
699 SolidEnergyLawParams solidEnergyLawParams_;
705 unsigned nTemperature_;
708 Scalar pressureLow_, pressureHigh_;
709 Scalar temperatureLow_, temperatureHigh_;