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>
48 #include <dune/grid/yaspgrid.hh>
49 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
51 #include <dune/common/version.hh>
52 #include <dune/common/fvector.hh>
53 #include <dune/common/fmatrix.hh>
61 template <
class TypeTag>
62 class Co2InjectionProblem;
64 namespace Co2Injection {
65 #include <opm/material/components/co2tables.inc>
70 namespace Opm::Properties {
77 template<
class TypeTag,
class MyTypeTag>
79 template<
class TypeTag,
class MyTypeTag>
81 template<
class TypeTag,
class MyTypeTag>
83 template<
class TypeTag,
class MyTypeTag>
85 template<
class TypeTag,
class MyTypeTag>
87 template<
class TypeTag,
class MyTypeTag>
90 template<
class TypeTag,
class MyTypeTag>
91 struct MaxDepth {
using type = UndefinedProperty; };
92 template<
class TypeTag,
class MyTypeTag>
94 template<
class TypeTag,
class MyTypeTag>
98 template<
class TypeTag>
99 struct Grid<TypeTag, TTag::Co2InjectionBaseProblem> {
using type = Dune::YaspGrid<2>; };
102 template<
class TypeTag>
103 struct Problem<TypeTag, TTag::Co2InjectionBaseProblem>
107 template<
class TypeTag>
108 struct FluidSystem<TypeTag, TTag::Co2InjectionBaseProblem>
111 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
112 using CO2Tables = Opm::Co2Injection::CO2Tables;
115 using type = Opm::BrineCO2FluidSystem<Scalar, CO2Tables>;
120 template<
class TypeTag>
121 struct MaterialLaw<TypeTag, TTag::Co2InjectionBaseProblem>
124 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
125 enum { liquidPhaseIdx = FluidSystem::liquidPhaseIdx };
126 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
128 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
129 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
130 FluidSystem::liquidPhaseIdx,
131 FluidSystem::gasPhaseIdx>;
135 using EffMaterialLaw = Opm::RegularizedBrooksCorey<Traits>;
139 using type = Opm::EffToAbsLaw<EffMaterialLaw>;
143 template<
class TypeTag>
144 struct ThermalConductionLaw<TypeTag, TTag::Co2InjectionBaseProblem>
147 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
148 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
152 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
156 template<
class TypeTag>
157 struct SolidEnergyLaw<TypeTag, TTag::Co2InjectionBaseProblem>
158 {
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
161 template<
class TypeTag>
162 struct LinearSolverSplice<TypeTag, TTag::Co2InjectionBaseProblem> {
using type = TTag::ParallelAmgLinearSolver; };
165 template<
class TypeTag>
166 struct NewtonWriteConvergence<TypeTag, TTag::Co2InjectionBaseProblem> {
static constexpr
bool value =
false; };
169 template<
class TypeTag>
170 struct EnableGravity<TypeTag, TTag::Co2InjectionBaseProblem> {
static constexpr
bool value =
true; };
173 template<
class TypeTag>
176 using type = GetPropType<TypeTag, Scalar>;
177 static constexpr type value = 3e7;
179 template<
class TypeTag>
182 using type = GetPropType<TypeTag, Scalar>;
183 static constexpr type value = 4e7;
185 template<
class TypeTag>
187 template<
class TypeTag>
190 using type = GetPropType<TypeTag, Scalar>;
191 static constexpr type value = 290;
193 template<
class TypeTag>
196 using type = GetPropType<TypeTag, Scalar>;
197 static constexpr type value = 500;
199 template<
class TypeTag>
202 template<
class TypeTag>
203 struct MaxDepth<TypeTag, TTag::Co2InjectionBaseProblem>
205 using type = GetPropType<TypeTag, Scalar>;
206 static constexpr type value = 2500;
208 template<
class TypeTag>
211 using type = GetPropType<TypeTag, Scalar>;
212 static constexpr type value = 293.15;
214 template<
class TypeTag>
215 struct SimulationName<TypeTag, TTag::Co2InjectionBaseProblem> {
static constexpr
auto value =
"co2injection"; };
218 template<
class TypeTag>
219 struct EndTime<TypeTag, TTag::Co2InjectionBaseProblem>
221 using type = GetPropType<TypeTag, Scalar>;
222 static constexpr type value = 1e4;
226 template<
class TypeTag>
227 struct InitialTimeStepSize<TypeTag, TTag::Co2InjectionBaseProblem>
229 using type = GetPropType<TypeTag, Scalar>;
230 static constexpr type value = 250;
234 template<
class TypeTag>
235 struct GridFile<TypeTag, TTag::Co2InjectionBaseProblem> {
static constexpr
auto value =
"data/co2injection.dgf"; };
262 template <
class TypeTag>
265 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
267 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
268 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
269 using GridView = GetPropType<TypeTag, Properties::GridView>;
270 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
272 enum { dim = GridView::dimension };
273 enum { dimWorld = GridView::dimensionworld };
276 using Indices = GetPropType<TypeTag, Properties::Indices>;
277 enum { numPhases = FluidSystem::numPhases };
278 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
279 enum { liquidPhaseIdx = FluidSystem::liquidPhaseIdx };
280 enum { CO2Idx = FluidSystem::CO2Idx };
281 enum { BrineIdx = FluidSystem::BrineIdx };
282 enum { conti0EqIdx = Indices::conti0EqIdx };
283 enum { contiCO2EqIdx = conti0EqIdx + CO2Idx };
285 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
286 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
287 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
288 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
289 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
290 using Model = GetPropType<TypeTag, Properties::Model>;
291 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
292 using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
293 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
294 using ThermalConductionLawParams =
typename ThermalConductionLaw::Params;
296 using Toolbox = Opm::MathToolbox<Evaluation>;
297 using CoordScalar =
typename GridView::ctype;
298 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
299 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
306 : ParentType(simulator)
314 ParentType::finishInit();
318 temperatureLow_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureLow);
319 temperatureHigh_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureHigh);
320 nTemperature_ = EWOMS_GET_PARAM(TypeTag,
unsigned, FluidSystemNumTemperature);
322 pressureLow_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemPressureLow);
323 pressureHigh_ = EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemPressureHigh);
324 nPressure_ = EWOMS_GET_PARAM(TypeTag,
unsigned, FluidSystemNumPressure);
326 maxDepth_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxDepth);
327 temperature_ = EWOMS_GET_PARAM(TypeTag, Scalar, Temperature);
331 FluidSystem::init(temperatureLow_,
338 fineLayerBottom_ = 22.0;
341 fineK_ = this->toDimMatrix_(1e-13);
342 coarseK_ = this->toDimMatrix_(1e-12);
346 coarsePorosity_ = 0.3;
349 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
350 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
351 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
352 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
355 fineMaterialParams_.setEntryPressure(1e4);
356 coarseMaterialParams_.setEntryPressure(5e3);
357 fineMaterialParams_.setLambda(2.0);
358 coarseMaterialParams_.setLambda(2.0);
360 fineMaterialParams_.finalize();
361 coarseMaterialParams_.finalize();
364 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
365 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
368 solidEnergyLawParams_.setSolidHeatCapacity(790.0
370 solidEnergyLawParams_.finalize();
378 ParentType::registerParameters();
380 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemTemperatureLow,
381 "The lower temperature [K] for tabulation of the "
383 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemTemperatureHigh,
384 "The upper temperature [K] for tabulation of the "
386 EWOMS_REGISTER_PARAM(TypeTag,
unsigned, FluidSystemNumTemperature,
387 "The number of intervals between the lower and "
388 "upper temperature");
390 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemPressureLow,
391 "The lower pressure [Pa] for tabulation of the "
393 EWOMS_REGISTER_PARAM(TypeTag, Scalar, FluidSystemPressureHigh,
394 "The upper pressure [Pa] for tabulation of the "
396 EWOMS_REGISTER_PARAM(TypeTag,
unsigned, FluidSystemNumPressure,
397 "The number of intervals between the lower and "
400 EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature,
401 "The temperature [K] in the reservoir");
402 EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxDepth,
403 "The maximum depth [m] of the reservoir");
404 EWOMS_REGISTER_PARAM(TypeTag, std::string, SimulationName,
405 "The name of the simulation used for the output "
419 std::ostringstream oss;
420 oss << EWOMS_GET_PARAM(TypeTag, std::string, SimulationName)
421 <<
"_" << Model::name();
422 if (getPropValue<TypeTag, Properties::EnableEnergy>())
424 oss <<
"_" << Model::discretizationName();
434 Scalar tol = this->model().newtonMethod().tolerance()*1e5;
435 this->model().checkConservativeness(tol);
438 PrimaryVariables storageL, storageG;
439 this->model().globalPhaseStorage(storageL, 0);
440 this->model().globalPhaseStorage(storageG, 1);
443 if (this->gridView().comm().rank() == 0) {
444 std::cout <<
"Storage: liquid=[" << storageL <<
"]"
445 <<
" gas=[" << storageG <<
"]\n" << std::flush;
453 template <
class Context>
454 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
456 const auto& pos = context.pos(spaceIdx, timeIdx);
457 if (inHighTemperatureRegion_(pos))
458 return temperature_ + 100;
465 template <
class Context>
467 unsigned timeIdx)
const
469 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
470 if (isFineMaterial_(pos))
478 template <
class Context>
479 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
481 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
482 if (isFineMaterial_(pos))
483 return finePorosity_;
484 return coarsePorosity_;
490 template <
class Context>
492 unsigned spaceIdx,
unsigned timeIdx)
const
494 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
495 if (isFineMaterial_(pos))
496 return fineMaterialParams_;
497 return coarseMaterialParams_;
505 template <
class Context>
506 const SolidEnergyLawParams&
510 {
return solidEnergyLawParams_; }
515 template <
class Context>
516 const ThermalConductionLawParams &
519 unsigned timeIdx)
const
521 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
522 if (isFineMaterial_(pos))
523 return fineThermalCondParams_;
524 return coarseThermalCondParams_;
537 template <
class Context>
538 void boundary(BoundaryRateVector& values,
const Context& context,
539 unsigned spaceIdx,
unsigned timeIdx)
const
541 const auto& pos = context.pos(spaceIdx, timeIdx);
542 if (onLeftBoundary_(pos)) {
543 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
544 initialFluidState_(fs, context, spaceIdx, timeIdx);
548 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
550 else if (onInlet_(pos)) {
551 RateVector massRate(0.0);
552 massRate[contiCO2EqIdx] = -1e-3;
554 using FluidState = Opm::ImmiscibleFluidState<Scalar, FluidSystem>;
556 fs.setSaturation(gasPhaseIdx, 1.0);
558 context.intensiveQuantities(spaceIdx, timeIdx).fluidState().pressure(gasPhaseIdx);
559 fs.setPressure(gasPhaseIdx, Toolbox::value(pg));
560 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
562 typename FluidSystem::template ParameterCache<Scalar> paramCache;
563 paramCache.updatePhase(fs, gasPhaseIdx);
564 Scalar h = FluidSystem::template enthalpy<FluidState, Scalar>(fs, paramCache, gasPhaseIdx);
567 values.setMassRate(massRate);
568 values.setEnthalpyRate(massRate[contiCO2EqIdx] * h);
585 template <
class Context>
586 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
587 unsigned timeIdx)
const
589 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
590 initialFluidState_(fs, context, spaceIdx, timeIdx);
595 values.assignNaive(fs);
604 template <
class Context>
609 { rate = Scalar(0.0); }
614 template <
class Context,
class Flu
idState>
615 void initialFluidState_(FluidState& fs,
616 const Context& context,
618 unsigned timeIdx)
const
620 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
625 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
630 fs.setSaturation(FluidSystem::liquidPhaseIdx, 1.0);
631 fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
636 Scalar densityL = FluidSystem::Brine::liquidDensity(temperature_, Scalar(1e5));
637 Scalar depth = maxDepth_ - pos[dim - 1];
638 Scalar pl = 1e5 - densityL * this->gravity()[dim - 1] * depth;
640 Scalar pC[numPhases];
642 MaterialLaw::capillaryPressures(pC, matParams, fs);
644 fs.setPressure(liquidPhaseIdx, pl + (pC[liquidPhaseIdx] - pC[liquidPhaseIdx]));
645 fs.setPressure(gasPhaseIdx, pl + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
650 fs.setMoleFraction(liquidPhaseIdx, CO2Idx, 0.005);
651 fs.setMoleFraction(liquidPhaseIdx, BrineIdx,
652 1.0 - fs.moleFraction(liquidPhaseIdx, CO2Idx));
654 typename FluidSystem::template ParameterCache<Scalar> paramCache;
655 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
656 CFRP::solve(fs, paramCache,
662 bool onLeftBoundary_(
const GlobalPosition& pos)
const
663 {
return pos[0] < eps_; }
665 bool onRightBoundary_(
const GlobalPosition& pos)
const
666 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
668 bool onInlet_(
const GlobalPosition& pos)
const
669 {
return onRightBoundary_(pos) && (5 < pos[1]) && (pos[1] < 15); }
671 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const
672 {
return (pos[0] > 20) && (pos[0] < 30) && (pos[1] > 5) && (pos[1] < 35); }
674 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
676 Scalar lambdaWater = 0.6;
677 Scalar lambdaGranite = 2.8;
679 Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
680 * std::pow(lambdaWater, poro);
681 Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
683 params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
684 params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
685 params.setVacuumLambda(lambdaDry);
688 bool isFineMaterial_(
const GlobalPosition& pos)
const
689 {
return pos[dim - 1] > fineLayerBottom_; }
693 Scalar fineLayerBottom_;
695 Scalar finePorosity_;
696 Scalar coarsePorosity_;
698 MaterialLawParams fineMaterialParams_;
699 MaterialLawParams coarseMaterialParams_;
701 ThermalConductionLawParams fineThermalCondParams_;
702 ThermalConductionLawParams coarseThermalCondParams_;
703 SolidEnergyLawParams solidEnergyLawParams_;
709 unsigned nTemperature_;
712 Scalar pressureLow_, pressureHigh_;
713 Scalar temperatureLow_, temperatureHigh_;
Problem where is injected under a low permeable layer at a depth of 2700m.
Definition: co2injectionproblem.hh:264
void finishInit()
Definition: co2injectionproblem.hh:312
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:479
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:538
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:586
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:491
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: co2injectionproblem.hh:605
std::string name() const
Definition: co2injectionproblem.hh:417
void endTimeStep()
Definition: co2injectionproblem.hh:431
const SolidEnergyLawParams & solidEnergyLawParams(const Context &, unsigned, unsigned) const
Return the parameters for the heat storage law of the rock.
Definition: co2injectionproblem.hh:507
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:517
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:454
Co2InjectionProblem(Simulator &simulator)
Definition: co2injectionproblem.hh:305
static void registerParameters()
Definition: co2injectionproblem.hh:376
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:466
Definition: co2injectionproblem.hh:82
Definition: co2injectionproblem.hh:88
Definition: co2injectionproblem.hh:80
Definition: co2injectionproblem.hh:78
Definition: co2injectionproblem.hh:86
Definition: co2injectionproblem.hh:84
Definition: co2injectionproblem.hh:91
Definition: co2injectionproblem.hh:95
Definition: co2injectionproblem.hh:73
Definition: co2injectionproblem.hh:93