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_ = Parameters::get<TypeTag, Properties::FluidSystemTemperatureLow>();
315 temperatureHigh_ = Parameters::get<TypeTag, Properties::FluidSystemTemperatureHigh>();
316 nTemperature_ = Parameters::get<TypeTag, Properties::FluidSystemNumTemperature>();
318 pressureLow_ = Parameters::get<TypeTag, Properties::FluidSystemPressureLow>();
319 pressureHigh_ = Parameters::get<TypeTag, Properties::FluidSystemPressureHigh>();
320 nPressure_ = Parameters::get<TypeTag, Properties::FluidSystemNumPressure>();
322 maxDepth_ = Parameters::get<TypeTag, Properties::MaxDepth>();
323 temperature_ = Parameters::get<TypeTag, Properties::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 Parameters::registerParam<TypeTag, Properties::FluidSystemTemperatureLow>
377 (
"The lower temperature [K] for tabulation of the fluid system");
378 Parameters::registerParam<TypeTag, Properties::FluidSystemTemperatureHigh>
379 (
"The upper temperature [K] for tabulation of the fluid system");
380 Parameters::registerParam<TypeTag, Properties::FluidSystemNumTemperature>
381 (
"The number of intervals between the lower and upper temperature");
382 Parameters::registerParam<TypeTag, Properties::FluidSystemPressureLow>
383 (
"The lower pressure [Pa] for tabulation of the fluid system");
384 Parameters::registerParam<TypeTag, Properties::FluidSystemPressureHigh>
385 (
"The upper pressure [Pa] for tabulation of the fluid system");
386 Parameters::registerParam<TypeTag, Properties::FluidSystemNumPressure>
387 (
"The number of intervals between the lower and upper pressure");
388 Parameters::registerParam<TypeTag, Properties::Temperature>
389 (
"The temperature [K] in the reservoir");
390 Parameters::registerParam<TypeTag, Properties::MaxDepth>
391 (
"The maximum depth [m] of the reservoir");
392 Parameters::registerParam<TypeTag, Properties::SimulationName>
393 (
"The name of the simulation used for the output files");
406 std::ostringstream oss;
407 oss << Parameters::get<TypeTag, Properties::SimulationName>()
408 <<
"_" << Model::name();
409 if (getPropValue<TypeTag, Properties::EnableEnergy>())
411 oss <<
"_" << Model::discretizationName();
421 Scalar tol = this->model().newtonMethod().tolerance()*1e5;
422 this->model().checkConservativeness(tol);
425 PrimaryVariables storageL, storageG;
426 this->model().globalPhaseStorage(storageL, 0);
427 this->model().globalPhaseStorage(storageG, 1);
430 if (this->gridView().comm().rank() == 0) {
431 std::cout <<
"Storage: liquid=[" << storageL <<
"]"
432 <<
" gas=[" << storageG <<
"]\n" << std::flush;
440 template <
class Context>
441 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
443 const auto& pos = context.pos(spaceIdx, timeIdx);
444 if (inHighTemperatureRegion_(pos))
445 return temperature_ + 100;
452 template <
class Context>
454 unsigned timeIdx)
const
456 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
457 if (isFineMaterial_(pos))
465 template <
class Context>
466 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
468 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
469 if (isFineMaterial_(pos))
470 return finePorosity_;
471 return coarsePorosity_;
477 template <
class Context>
479 unsigned spaceIdx,
unsigned timeIdx)
const
481 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
482 if (isFineMaterial_(pos))
483 return fineMaterialParams_;
484 return coarseMaterialParams_;
492 template <
class Context>
493 const SolidEnergyLawParams&
497 {
return solidEnergyLawParams_; }
502 template <
class Context>
503 const ThermalConductionLawParams &
506 unsigned timeIdx)
const
508 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
509 if (isFineMaterial_(pos))
510 return fineThermalCondParams_;
511 return coarseThermalCondParams_;
524 template <
class Context>
525 void boundary(BoundaryRateVector& values,
const Context& context,
526 unsigned spaceIdx,
unsigned timeIdx)
const
528 const auto& pos = context.pos(spaceIdx, timeIdx);
529 if (onLeftBoundary_(pos)) {
530 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
531 initialFluidState_(fs, context, spaceIdx, timeIdx);
535 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
537 else if (onInlet_(pos)) {
538 RateVector massRate(0.0);
539 massRate[contiCO2EqIdx] = -1e-3;
541 using FluidState = Opm::ImmiscibleFluidState<Scalar, FluidSystem>;
543 fs.setSaturation(gasPhaseIdx, 1.0);
545 context.intensiveQuantities(spaceIdx, timeIdx).fluidState().pressure(gasPhaseIdx);
546 fs.setPressure(gasPhaseIdx, Toolbox::value(pg));
547 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
549 typename FluidSystem::template ParameterCache<Scalar> paramCache;
550 paramCache.updatePhase(fs, gasPhaseIdx);
551 Scalar h = FluidSystem::template enthalpy<FluidState, Scalar>(fs, paramCache, gasPhaseIdx);
554 values.setMassRate(massRate);
555 values.setEnthalpyRate(massRate[contiCO2EqIdx] * h);
572 template <
class Context>
573 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
574 unsigned timeIdx)
const
576 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
577 initialFluidState_(fs, context, spaceIdx, timeIdx);
582 values.assignNaive(fs);
591 template <
class Context>
596 { rate = Scalar(0.0); }
601 template <
class Context,
class Flu
idState>
602 void initialFluidState_(FluidState& fs,
603 const Context& context,
605 unsigned timeIdx)
const
607 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
612 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
617 fs.setSaturation(FluidSystem::liquidPhaseIdx, 1.0);
618 fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
623 Scalar densityL = FluidSystem::Brine::liquidDensity(temperature_, Scalar(1e5));
624 Scalar depth = maxDepth_ - pos[dim - 1];
625 Scalar pl = 1e5 - densityL * this->gravity()[dim - 1] * depth;
627 Scalar pC[numPhases];
629 MaterialLaw::capillaryPressures(pC, matParams, fs);
631 fs.setPressure(liquidPhaseIdx, pl + (pC[liquidPhaseIdx] - pC[liquidPhaseIdx]));
632 fs.setPressure(gasPhaseIdx, pl + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
637 fs.setMoleFraction(liquidPhaseIdx, CO2Idx, 0.005);
638 fs.setMoleFraction(liquidPhaseIdx, BrineIdx,
639 1.0 - fs.moleFraction(liquidPhaseIdx, CO2Idx));
641 typename FluidSystem::template ParameterCache<Scalar> paramCache;
642 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
643 CFRP::solve(fs, paramCache,
649 bool onLeftBoundary_(
const GlobalPosition& pos)
const
650 {
return pos[0] < eps_; }
652 bool onRightBoundary_(
const GlobalPosition& pos)
const
653 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
655 bool onInlet_(
const GlobalPosition& pos)
const
656 {
return onRightBoundary_(pos) && (5 < pos[1]) && (pos[1] < 15); }
658 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const
659 {
return (pos[0] > 20) && (pos[0] < 30) && (pos[1] > 5) && (pos[1] < 35); }
661 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
663 Scalar lambdaWater = 0.6;
664 Scalar lambdaGranite = 2.8;
666 Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
667 * std::pow(lambdaWater, poro);
668 Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
670 params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
671 params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
672 params.setVacuumLambda(lambdaDry);
675 bool isFineMaterial_(
const GlobalPosition& pos)
const
676 {
return pos[dim - 1] > fineLayerBottom_; }
680 Scalar fineLayerBottom_;
682 Scalar finePorosity_;
683 Scalar coarsePorosity_;
685 MaterialLawParams fineMaterialParams_;
686 MaterialLawParams coarseMaterialParams_;
688 ThermalConductionLawParams fineThermalCondParams_;
689 ThermalConductionLawParams coarseThermalCondParams_;
690 SolidEnergyLawParams solidEnergyLawParams_;
696 unsigned nTemperature_;
699 Scalar pressureLow_, pressureHigh_;
700 Scalar temperatureLow_, temperatureHigh_;