28 #ifndef EWOMS_CO2_INJECTION_PROBLEM_HH
29 #define EWOMS_CO2_INJECTION_PROBLEM_HH
31 #include <opm/models/immiscible/immisciblemodel.hh>
32 #include <opm/simulators/linalg/parallelamgbackend.hh>
34 #include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
35 #include <opm/material/fluidsystems/BrineCO2FluidSystem.hpp>
36 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
37 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
38 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
39 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
40 #include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
41 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
42 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
43 #include <opm/material/thermal/SomertonThermalConductionLaw.hpp>
44 #include <opm/material/thermal/ConstantSolidHeatCapLaw.hpp>
45 #include <opm/material/binarycoefficients/Brine_CO2.hpp>
46 #include <opm/material/common/UniformTabulated2DFunction.hpp>
47 #include <opm/material/common/Unused.hpp>
49 #include <dune/grid/yaspgrid.hh>
50 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
52 #include <dune/common/version.hh>
53 #include <dune/common/fvector.hh>
54 #include <dune/common/fmatrix.hh>
62 template <
class TypeTag>
63 class Co2InjectionProblem;
65 namespace Co2Injection {
66 #include <opm/material/components/co2tables.inc>
71 namespace Opm::Properties {
78 template<
class TypeTag,
class MyTypeTag>
80 template<
class TypeTag,
class MyTypeTag>
82 template<
class TypeTag,
class MyTypeTag>
84 template<
class TypeTag,
class MyTypeTag>
86 template<
class TypeTag,
class MyTypeTag>
88 template<
class TypeTag,
class MyTypeTag>
91 template<
class TypeTag,
class MyTypeTag>
92 struct MaxDepth {
using type = UndefinedProperty; };
93 template<
class TypeTag,
class MyTypeTag>
95 template<
class TypeTag,
class MyTypeTag>
99 template<
class TypeTag>
100 struct Grid<TypeTag, TTag::Co2InjectionBaseProblem> {
using type = Dune::YaspGrid<2>; };
103 template<
class TypeTag>
104 struct Problem<TypeTag, TTag::Co2InjectionBaseProblem>
108 template<
class TypeTag>
109 struct FluidSystem<TypeTag, TTag::Co2InjectionBaseProblem>
112 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
113 using CO2Tables = Opm::Co2Injection::CO2Tables;
116 using type = Opm::BrineCO2FluidSystem<Scalar, CO2Tables>;
121 template<
class TypeTag>
122 struct MaterialLaw<TypeTag, TTag::Co2InjectionBaseProblem>
125 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
126 enum { liquidPhaseIdx = FluidSystem::liquidPhaseIdx };
127 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
129 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
130 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
131 FluidSystem::liquidPhaseIdx,
132 FluidSystem::gasPhaseIdx>;
136 using EffMaterialLaw = Opm::RegularizedBrooksCorey<Traits>;
140 using type = Opm::EffToAbsLaw<EffMaterialLaw>;
144 template<
class TypeTag>
145 struct ThermalConductionLaw<TypeTag, TTag::Co2InjectionBaseProblem>
148 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
149 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
153 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
157 template<
class TypeTag>
158 struct SolidEnergyLaw<TypeTag, TTag::Co2InjectionBaseProblem>
159 {
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
162 template<
class TypeTag>
163 struct LinearSolverSplice<TypeTag, TTag::Co2InjectionBaseProblem> {
using type = TTag::ParallelAmgLinearSolver; };
166 template<
class TypeTag>
167 struct NewtonWriteConvergence<TypeTag, TTag::Co2InjectionBaseProblem> {
static constexpr
bool value =
false; };
170 template<
class TypeTag>
171 struct EnableGravity<TypeTag, TTag::Co2InjectionBaseProblem> {
static constexpr
bool value =
true; };
174 template<
class TypeTag>
177 using type = GetPropType<TypeTag, Scalar>;
178 static constexpr type value = 3e7;
180 template<
class TypeTag>
183 using type = GetPropType<TypeTag, Scalar>;
184 static constexpr type value = 4e7;
186 template<
class TypeTag>
188 template<
class TypeTag>
191 using type = GetPropType<TypeTag, Scalar>;
192 static constexpr type value = 290;
194 template<
class TypeTag>
197 using type = GetPropType<TypeTag, Scalar>;
198 static constexpr type value = 500;
200 template<
class TypeTag>
203 template<
class TypeTag>
204 struct MaxDepth<TypeTag, TTag::Co2InjectionBaseProblem>
206 using type = GetPropType<TypeTag, Scalar>;
207 static constexpr type value = 2500;
209 template<
class TypeTag>
212 using type = GetPropType<TypeTag, Scalar>;
213 static constexpr type value = 293.15;
215 template<
class TypeTag>
216 struct SimulationName<TypeTag, TTag::Co2InjectionBaseProblem> {
static constexpr
auto value =
"co2injection"; };
219 template<
class TypeTag>
220 struct EndTime<TypeTag, TTag::Co2InjectionBaseProblem>
222 using type = GetPropType<TypeTag, Scalar>;
223 static constexpr type value = 1e4;
227 template<
class TypeTag>
228 struct InitialTimeStepSize<TypeTag, TTag::Co2InjectionBaseProblem>
230 using type = GetPropType<TypeTag, Scalar>;
231 static constexpr type value = 250;
235 template<
class TypeTag>
236 struct GridFile<TypeTag, TTag::Co2InjectionBaseProblem> {
static constexpr
auto value =
"data/co2injection.dgf"; };
263 template <
class TypeTag>
266 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
268 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
269 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
270 using GridView = GetPropType<TypeTag, Properties::GridView>;
271 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
273 enum { dim = GridView::dimension };
274 enum { dimWorld = GridView::dimensionworld };
277 using Indices = GetPropType<TypeTag, Properties::Indices>;
278 enum { numPhases = FluidSystem::numPhases };
279 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
280 enum { liquidPhaseIdx = FluidSystem::liquidPhaseIdx };
281 enum { CO2Idx = FluidSystem::CO2Idx };
282 enum { BrineIdx = FluidSystem::BrineIdx };
283 enum { conti0EqIdx = Indices::conti0EqIdx };
284 enum { contiCO2EqIdx = conti0EqIdx + CO2Idx };
286 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
287 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
288 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
289 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
290 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
291 using Model = GetPropType<TypeTag, Properties::Model>;
292 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
293 using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
294 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
295 using ThermalConductionLawParams =
typename ThermalConductionLaw::Params;
297 using Toolbox = Opm::MathToolbox<Evaluation>;
298 using CoordScalar =
typename GridView::ctype;
299 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
300 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
307 : ParentType(simulator)
315 ParentType::finishInit();
319 temperatureLow_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureLow);
320 temperatureHigh_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureHigh);
321 nTemperature_ = EWOMS_GET_PARAM(TypeTag,
unsigned, FluidSystemNumTemperature);
323 pressureLow_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemPressureLow);
324 pressureHigh_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemPressureHigh);
325 nPressure_ = EWOMS_GET_PARAM(TypeTag,
unsigned, FluidSystemNumPressure);
327 maxDepth_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxDepth);
328 temperature_ = EWOMS_GET_PARAM(TypeTag, Scalar, Temperature);
332 FluidSystem::init(temperatureLow_,
339 fineLayerBottom_ = 22.0;
342 fineK_ = this->toDimMatrix_(1e-13);
343 coarseK_ = this->toDimMatrix_(1e-12);
347 coarsePorosity_ = 0.3;
350 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
351 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
352 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
353 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
356 fineMaterialParams_.setEntryPressure(1e4);
357 coarseMaterialParams_.setEntryPressure(5e3);
358 fineMaterialParams_.setLambda(2.0);
359 coarseMaterialParams_.setLambda(2.0);
361 fineMaterialParams_.finalize();
362 coarseMaterialParams_.finalize();
365 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
366 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
369 solidEnergyLawParams_.setSolidHeatCapacity(790.0
371 solidEnergyLawParams_.finalize();
379 ParentType::registerParameters();
381 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemTemperatureLow,
382 "The lower temperature [K] for tabulation of the "
384 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemTemperatureHigh,
385 "The upper temperature [K] for tabulation of the "
387 EWOMS_REGISTER_PARAM(TypeTag,
unsigned, FluidSystemNumTemperature,
388 "The number of intervals between the lower and "
389 "upper temperature");
391 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemPressureLow,
392 "The lower pressure [Pa] for tabulation of the "
394 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemPressureHigh,
395 "The upper pressure [Pa] for tabulation of the "
397 EWOMS_REGISTER_PARAM(TypeTag,
unsigned, FluidSystemNumPressure,
398 "The number of intervals between the lower and "
401 EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature,
402 "The temperature [K] in the reservoir");
403 EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxDepth,
404 "The maximum depth [m] of the reservoir");
405 EWOMS_REGISTER_PARAM(TypeTag, std::string, SimulationName,
406 "The name of the simulation used for the output "
420 std::ostringstream oss;
421 oss << EWOMS_GET_PARAM(TypeTag, std::string, SimulationName)
422 <<
"_" << Model::name();
423 if (getPropValue<TypeTag, Properties::EnableEnergy>())
425 oss <<
"_" << Model::discretizationName();
435 Scalar tol = this->model().newtonMethod().tolerance()*1e5;
436 this->model().checkConservativeness(tol);
439 PrimaryVariables storageL, storageG;
440 this->model().globalPhaseStorage(storageL, 0);
441 this->model().globalPhaseStorage(storageG, 1);
444 if (this->gridView().comm().rank() == 0) {
445 std::cout <<
"Storage: liquid=[" << storageL <<
"]"
446 <<
" gas=[" << storageG <<
"]\n" << std::flush;
454 template <
class Context>
455 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
457 const auto& pos = context.pos(spaceIdx, timeIdx);
458 if (inHighTemperatureRegion_(pos))
459 return temperature_ + 100;
466 template <
class Context>
468 unsigned timeIdx)
const
470 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
471 if (isFineMaterial_(pos))
479 template <
class Context>
480 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
482 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
483 if (isFineMaterial_(pos))
484 return finePorosity_;
485 return coarsePorosity_;
491 template <
class Context>
493 unsigned spaceIdx,
unsigned timeIdx)
const
495 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
496 if (isFineMaterial_(pos))
497 return fineMaterialParams_;
498 return coarseMaterialParams_;
506 template <
class Context>
507 const SolidEnergyLawParams&
509 unsigned spaceIdx OPM_UNUSED,
510 unsigned timeIdx OPM_UNUSED)
const
511 {
return solidEnergyLawParams_; }
516 template <
class Context>
517 const ThermalConductionLawParams &
520 unsigned timeIdx)
const
522 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
523 if (isFineMaterial_(pos))
524 return fineThermalCondParams_;
525 return coarseThermalCondParams_;
538 template <
class Context>
539 void boundary(BoundaryRateVector& values,
const Context& context,
540 unsigned spaceIdx,
unsigned timeIdx)
const
542 const auto& pos = context.pos(spaceIdx, timeIdx);
543 if (onLeftBoundary_(pos)) {
544 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
545 initialFluidState_(fs, context, spaceIdx, timeIdx);
549 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
551 else if (onInlet_(pos)) {
552 RateVector massRate(0.0);
553 massRate[contiCO2EqIdx] = -1e-3;
555 using FluidState = Opm::ImmiscibleFluidState<Scalar, FluidSystem>;
557 fs.setSaturation(gasPhaseIdx, 1.0);
559 context.intensiveQuantities(spaceIdx, timeIdx).fluidState().pressure(gasPhaseIdx);
560 fs.setPressure(gasPhaseIdx, Toolbox::value(pg));
561 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
563 typename FluidSystem::template ParameterCache<Scalar> paramCache;
564 paramCache.updatePhase(fs, gasPhaseIdx);
565 Scalar h = FluidSystem::template enthalpy<FluidState, Scalar>(fs, paramCache, gasPhaseIdx);
568 values.setMassRate(massRate);
569 values.setEnthalpyRate(massRate[contiCO2EqIdx] * h);
586 template <
class Context>
587 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
588 unsigned timeIdx)
const
590 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
591 initialFluidState_(fs, context, spaceIdx, timeIdx);
596 values.assignNaive(fs);
605 template <
class Context>
607 const Context& context OPM_UNUSED,
608 unsigned spaceIdx OPM_UNUSED,
609 unsigned timeIdx OPM_UNUSED)
const
610 { rate = Scalar(0.0); }
615 template <
class Context,
class Flu
idState>
616 void initialFluidState_(FluidState& fs,
617 const Context& context,
619 unsigned timeIdx)
const
621 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
626 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
631 fs.setSaturation(FluidSystem::liquidPhaseIdx, 1.0);
632 fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
637 Scalar densityL = FluidSystem::Brine::liquidDensity(temperature_, Scalar(1e5));
638 Scalar depth = maxDepth_ - pos[dim - 1];
639 Scalar pl = 1e5 - densityL * this->gravity()[dim - 1] * depth;
641 Scalar pC[numPhases];
643 MaterialLaw::capillaryPressures(pC, matParams, fs);
645 fs.setPressure(liquidPhaseIdx, pl + (pC[liquidPhaseIdx] - pC[liquidPhaseIdx]));
646 fs.setPressure(gasPhaseIdx, pl + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
651 fs.setMoleFraction(liquidPhaseIdx, CO2Idx, 0.005);
652 fs.setMoleFraction(liquidPhaseIdx, BrineIdx,
653 1.0 - fs.moleFraction(liquidPhaseIdx, CO2Idx));
655 typename FluidSystem::template ParameterCache<Scalar> paramCache;
656 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
657 CFRP::solve(fs, paramCache,
663 bool onLeftBoundary_(
const GlobalPosition& pos)
const
664 {
return pos[0] < eps_; }
666 bool onRightBoundary_(
const GlobalPosition& pos)
const
667 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
669 bool onInlet_(
const GlobalPosition& pos)
const
670 {
return onRightBoundary_(pos) && (5 < pos[1]) && (pos[1] < 15); }
672 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const
673 {
return (pos[0] > 20) && (pos[0] < 30) && (pos[1] > 5) && (pos[1] < 35); }
675 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
677 Scalar lambdaWater = 0.6;
678 Scalar lambdaGranite = 2.8;
680 Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
681 * std::pow(lambdaWater, poro);
682 Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
684 params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
685 params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
686 params.setVacuumLambda(lambdaDry);
689 bool isFineMaterial_(
const GlobalPosition& pos)
const
690 {
return pos[dim - 1] > fineLayerBottom_; }
694 Scalar fineLayerBottom_;
696 Scalar finePorosity_;
697 Scalar coarsePorosity_;
699 MaterialLawParams fineMaterialParams_;
700 MaterialLawParams coarseMaterialParams_;
702 ThermalConductionLawParams fineThermalCondParams_;
703 ThermalConductionLawParams coarseThermalCondParams_;
704 SolidEnergyLawParams solidEnergyLawParams_;
710 unsigned nTemperature_;
713 Scalar pressureLow_, pressureHigh_;
714 Scalar temperatureLow_, temperatureHigh_;
Problem where is injected under a low permeable layer at a depth of 2700m.
Definition: co2injectionproblem.hh:265
void finishInit()
Definition: co2injectionproblem.hh:313
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:480
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:539
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:587
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:492
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: co2injectionproblem.hh:606
std::string name() const
Definition: co2injectionproblem.hh:418
void endTimeStep()
Definition: co2injectionproblem.hh:432
const SolidEnergyLawParams & solidEnergyLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Return the parameters for the heat storage law of the rock.
Definition: co2injectionproblem.hh:508
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:518
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:455
Co2InjectionProblem(Simulator &simulator)
Definition: co2injectionproblem.hh:306
static void registerParameters()
Definition: co2injectionproblem.hh:377
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:467
Definition: co2injectionproblem.hh:83
Definition: co2injectionproblem.hh:89
Definition: co2injectionproblem.hh:81
Definition: co2injectionproblem.hh:79
Definition: co2injectionproblem.hh:87
Definition: co2injectionproblem.hh:85
Definition: co2injectionproblem.hh:92
Definition: co2injectionproblem.hh:96
Definition: co2injectionproblem.hh:74
Definition: co2injectionproblem.hh:94