28#ifndef EWOMS_CO2_INJECTION_PROBLEM_HH
29#define EWOMS_CO2_INJECTION_PROBLEM_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>
62template <
class TypeTag>
63class Co2InjectionProblem;
75template<
class TypeTag>
76struct Grid<TypeTag, TTag::Co2InjectionBaseProblem> {
using type = Dune::YaspGrid<2>; };
79template<
class TypeTag>
80struct Problem<TypeTag, TTag::Co2InjectionBaseProblem>
84template<
class TypeTag>
91 using type = Opm::BrineCO2FluidSystem<Scalar>;
96template<
class TypeTag>
101 enum { liquidPhaseIdx = FluidSystem::liquidPhaseIdx };
102 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
105 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
106 FluidSystem::liquidPhaseIdx,
107 FluidSystem::gasPhaseIdx>;
111 using EffMaterialLaw = Opm::RegularizedBrooksCorey<Traits>;
115 using type = Opm::EffToAbsLaw<EffMaterialLaw>;
119template<
class TypeTag>
128 using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
132template<
class TypeTag>
134{
using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
137template<
class TypeTag>
147template<
class Scalar>
150template<
class Scalar>
153template<
class Scalar>
156template<
class Scalar>
159template<
class Scalar>
164template<
class Scalar>
195template <
class TypeTag>
205 enum { dim = GridView::dimension };
206 enum { dimWorld = GridView::dimensionworld };
210 enum { numPhases = FluidSystem::numPhases };
211 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
212 enum { liquidPhaseIdx = FluidSystem::liquidPhaseIdx };
213 enum { CO2Idx = FluidSystem::CO2Idx };
214 enum { BrineIdx = FluidSystem::BrineIdx };
215 enum { conti0EqIdx = Indices::conti0EqIdx };
216 enum { contiCO2EqIdx = conti0EqIdx + CO2Idx };
227 using ThermalConductionLawParams =
typename ThermalConductionLaw::Params;
229 using Toolbox = Opm::MathToolbox<Evaluation>;
230 using CoordScalar =
typename GridView::ctype;
231 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
232 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
239 : ParentType(simulator)
247 ParentType::finishInit();
251 temperatureLow_ = Parameters::Get<Parameters::FluidSystemTemperatureLow<Scalar>>();
252 temperatureHigh_ = Parameters::Get<Parameters::FluidSystemTemperatureHigh<Scalar>>();
253 nTemperature_ = Parameters::Get<Parameters::FluidSystemNumTemperature>();
255 pressureLow_ = Parameters::Get<Parameters::FluidSystemPressureLow<Scalar>>();
256 pressureHigh_ = Parameters::Get<Parameters::FluidSystemPressureHigh<Scalar>>();
257 nPressure_ = Parameters::Get<Parameters::FluidSystemNumPressure>();
259 maxDepth_ = Parameters::Get<Parameters::MaxDepth<Scalar>>();
260 temperature_ = Parameters::Get<Parameters::Temperature<Scalar>>();
264 FluidSystem::init(temperatureLow_,
271 fineLayerBottom_ = 22.0;
274 fineK_ = this->toDimMatrix_(1e-13);
275 coarseK_ = this->toDimMatrix_(1e-12);
279 coarsePorosity_ = 0.3;
282 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
283 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
284 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
285 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
288 fineMaterialParams_.setEntryPressure(1e4);
289 coarseMaterialParams_.setEntryPressure(5e3);
290 fineMaterialParams_.setLambda(2.0);
291 coarseMaterialParams_.setLambda(2.0);
293 fineMaterialParams_.finalize();
294 coarseMaterialParams_.finalize();
297 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
298 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
301 solidEnergyLawParams_.setSolidHeatCapacity(790.0
303 solidEnergyLawParams_.finalize();
311 ParentType::registerParameters();
313 Parameters::Register<Parameters::FluidSystemTemperatureLow<Scalar>>
314 (
"The lower temperature [K] for tabulation of the fluid system");
315 Parameters::Register<Parameters::FluidSystemTemperatureHigh<Scalar>>
316 (
"The upper temperature [K] for tabulation of the fluid system");
317 Parameters::Register<Parameters::FluidSystemNumTemperature>
318 (
"The number of intervals between the lower and upper temperature");
319 Parameters::Register<Parameters::FluidSystemPressureLow<Scalar>>
320 (
"The lower pressure [Pa] for tabulation of the fluid system");
321 Parameters::Register<Parameters::FluidSystemPressureHigh<Scalar>>
322 (
"The upper pressure [Pa] for tabulation of the fluid system");
323 Parameters::Register<Parameters::FluidSystemNumPressure>
324 (
"The number of intervals between the lower and upper pressure");
325 Parameters::Register<Parameters::Temperature<Scalar>>
326 (
"The temperature [K] in the reservoir");
327 Parameters::Register<Parameters::MaxDepth<Scalar>>
328 (
"The maximum depth [m] of the reservoir");
329 Parameters::Register<Parameters::SimulationName>
330 (
"The name of the simulation used for the output files");
332 Parameters::SetDefault<Parameters::GridFile>(
"data/co2injection.dgf");
333 Parameters::SetDefault<Parameters::EndTime<Scalar>>(1e4);
334 Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(250);
336 Parameters::SetDefault<Parameters::EnableGravity>(
true);
349 std::ostringstream oss;
350 oss << Parameters::Get<Parameters::SimulationName>()
351 <<
"_" << Model::name();
352 if (getPropValue<TypeTag, Properties::EnableEnergy>())
354 oss <<
"_" << Model::discretizationName();
364 Scalar tol = this->model().newtonMethod().tolerance()*1e5;
365 this->model().checkConservativeness(tol);
368 PrimaryVariables storageL, storageG;
369 this->model().globalPhaseStorage(storageL, 0);
370 this->model().globalPhaseStorage(storageG, 1);
373 if (this->gridView().comm().rank() == 0) {
374 std::cout <<
"Storage: liquid=[" << storageL <<
"]"
375 <<
" gas=[" << storageG <<
"]\n" << std::flush;
383 template <
class Context>
384 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
386 const auto& pos = context.pos(spaceIdx, timeIdx);
387 if (inHighTemperatureRegion_(pos))
388 return temperature_ + 100;
395 template <
class Context>
397 unsigned timeIdx)
const
399 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
400 if (isFineMaterial_(pos))
408 template <
class Context>
409 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
411 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
412 if (isFineMaterial_(pos))
413 return finePorosity_;
414 return coarsePorosity_;
420 template <
class Context>
422 unsigned spaceIdx,
unsigned timeIdx)
const
424 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
425 if (isFineMaterial_(pos))
426 return fineMaterialParams_;
427 return coarseMaterialParams_;
435 template <
class Context>
436 const SolidEnergyLawParams&
440 {
return solidEnergyLawParams_; }
445 template <
class Context>
446 const ThermalConductionLawParams &
449 unsigned timeIdx)
const
451 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
452 if (isFineMaterial_(pos))
453 return fineThermalCondParams_;
454 return coarseThermalCondParams_;
467 template <
class Context>
468 void boundary(BoundaryRateVector& values,
const Context& context,
469 unsigned spaceIdx,
unsigned timeIdx)
const
471 const auto& pos = context.pos(spaceIdx, timeIdx);
472 if (onLeftBoundary_(pos)) {
473 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
474 initialFluidState_(fs, context, spaceIdx, timeIdx);
478 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
480 else if (onInlet_(pos)) {
481 RateVector massRate(0.0);
482 massRate[contiCO2EqIdx] = -1e-3;
484 using FluidState = Opm::ImmiscibleFluidState<Scalar, FluidSystem>;
486 fs.setSaturation(gasPhaseIdx, 1.0);
488 context.intensiveQuantities(spaceIdx, timeIdx).fluidState().pressure(gasPhaseIdx);
489 fs.setPressure(gasPhaseIdx, Toolbox::value(pg));
490 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
492 typename FluidSystem::template ParameterCache<Scalar> paramCache;
493 paramCache.updatePhase(fs, gasPhaseIdx);
494 Scalar h = FluidSystem::template enthalpy<FluidState, Scalar>(fs, paramCache, gasPhaseIdx);
497 values.setMassRate(massRate);
498 values.setEnthalpyRate(massRate[contiCO2EqIdx] * h);
515 template <
class Context>
516 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
517 unsigned timeIdx)
const
519 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
520 initialFluidState_(fs, context, spaceIdx, timeIdx);
525 values.assignNaive(fs);
534 template <
class Context>
539 { rate = Scalar(0.0); }
544 template <
class Context,
class Flu
idState>
545 void initialFluidState_(FluidState& fs,
546 const Context& context,
548 unsigned timeIdx)
const
550 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
555 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
560 fs.setSaturation(FluidSystem::liquidPhaseIdx, 1.0);
561 fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
566 Scalar densityL = FluidSystem::Brine::liquidDensity(temperature_, Scalar(1e5));
567 Scalar depth = maxDepth_ - pos[dim - 1];
568 Scalar pl = 1e5 - densityL * this->gravity()[dim - 1] * depth;
570 Scalar pC[numPhases];
572 MaterialLaw::capillaryPressures(pC, matParams, fs);
574 fs.setPressure(liquidPhaseIdx, pl + (pC[liquidPhaseIdx] - pC[liquidPhaseIdx]));
575 fs.setPressure(gasPhaseIdx, pl + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
580 fs.setMoleFraction(liquidPhaseIdx, CO2Idx, 0.005);
581 fs.setMoleFraction(liquidPhaseIdx, BrineIdx,
582 1.0 - fs.moleFraction(liquidPhaseIdx, CO2Idx));
584 typename FluidSystem::template ParameterCache<Scalar> paramCache;
585 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
586 CFRP::solve(fs, paramCache,
592 bool onLeftBoundary_(
const GlobalPosition& pos)
const
593 {
return pos[0] < eps_; }
595 bool onRightBoundary_(
const GlobalPosition& pos)
const
596 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
598 bool onInlet_(
const GlobalPosition& pos)
const
599 {
return onRightBoundary_(pos) && (5 < pos[1]) && (pos[1] < 15); }
601 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const
602 {
return (pos[0] > 20) && (pos[0] < 30) && (pos[1] > 5) && (pos[1] < 35); }
604 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
606 Scalar lambdaWater = 0.6;
607 Scalar lambdaGranite = 2.8;
609 Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
610 * std::pow(lambdaWater, poro);
611 Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
613 params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
614 params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
615 params.setVacuumLambda(lambdaDry);
618 bool isFineMaterial_(
const GlobalPosition& pos)
const
619 {
return pos[dim - 1] > fineLayerBottom_; }
623 Scalar fineLayerBottom_;
625 Scalar finePorosity_;
626 Scalar coarsePorosity_;
628 MaterialLawParams fineMaterialParams_;
629 MaterialLawParams coarseMaterialParams_;
631 ThermalConductionLawParams fineThermalCondParams_;
632 ThermalConductionLawParams coarseThermalCondParams_;
633 SolidEnergyLawParams solidEnergyLawParams_;
639 unsigned nTemperature_;
642 Scalar pressureLow_, pressureHigh_;
643 Scalar temperatureLow_, temperatureHigh_;
Problem where is injected under a low permeable layer at a depth of 2700m.
Definition: co2injectionproblem.hh:197
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: co2injectionproblem.hh:245
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:409
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:447
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:396
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: co2injectionproblem.hh:468
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: co2injectionproblem.hh:516
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:421
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: co2injectionproblem.hh:535
const SolidEnergyLawParams & solidEnergyLawParams(const Context &, unsigned, unsigned) const
Return the parameters for the heat storage law of the rock.
Definition: co2injectionproblem.hh:437
std::string name() const
The problem name.
Definition: co2injectionproblem.hh:347
void endTimeStep()
Called by the simulator after each time integration.
Definition: co2injectionproblem.hh:361
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:384
Co2InjectionProblem(Simulator &simulator)
Definition: co2injectionproblem.hh:238
static void registerParameters()
Definition: co2injectionproblem.hh:309
Definition: blackoilnewtonmethodparams.hpp:31
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:235
double Co2InjectionTolerance
Definition: co2injectionproblem.hh:171
Definition: co2injectionproblem.hh:144
static constexpr unsigned value
Definition: co2injectionproblem.hh:144
Definition: co2injectionproblem.hh:145
static constexpr unsigned value
Definition: co2injectionproblem.hh:145
Definition: co2injectionproblem.hh:148
static constexpr Scalar value
Definition: co2injectionproblem.hh:148
Definition: co2injectionproblem.hh:151
static constexpr Scalar value
Definition: co2injectionproblem.hh:151
Definition: co2injectionproblem.hh:154
static constexpr Scalar value
Definition: co2injectionproblem.hh:154
Definition: co2injectionproblem.hh:157
static constexpr Scalar value
Definition: co2injectionproblem.hh:157
Definition: co2injectionproblem.hh:160
static constexpr Scalar value
Definition: co2injectionproblem.hh:160
Definition: co2injectionproblem.hh:162
static constexpr auto value
Definition: co2injectionproblem.hh:162
Definition: co2injectionproblem.hh:165
static constexpr Scalar value
Definition: co2injectionproblem.hh:165
Opm::BrineCO2FluidSystem< Scalar > type
Definition: co2injectionproblem.hh:91
The fluid systems including the information about the phases.
Definition: multiphasebaseproperties.hh:69
Dune::YaspGrid< 2 > type
Definition: co2injectionproblem.hh:76
The type of the DUNE grid.
Definition: basicproperties.hh:100
Definition: fvbaseproperties.hh:53
Opm::EffToAbsLaw< EffMaterialLaw > type
Definition: co2injectionproblem.hh:115
The material law which ought to be used (extracted from the spatial parameters)
Definition: multiphasebaseproperties.hh:51
The type of the problem.
Definition: fvbaseproperties.hh:81
Opm::ConstantSolidHeatCapLaw< GetPropType< TypeTag, Properties::Scalar > > type
Definition: co2injectionproblem.hh:134
The material law for the energy stored in the solid matrix.
Definition: multiphasebaseproperties.hh:57
Definition: co2injectionproblem.hh:71
Definition: parallelamgbackend.hh:58
Opm::SomertonThermalConductionLaw< FluidSystem, Scalar > type
Definition: co2injectionproblem.hh:128
The material law for thermal conduction.
Definition: multiphasebaseproperties.hh:63