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/heatconduction/Somerton.hpp> 44 #include <opm/material/binarycoefficients/Brine_CO2.hpp> 45 #include <opm/material/common/UniformTabulated2DFunction.hpp> 46 #include <opm/common/Unused.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> 69 namespace Properties {
85 SET_TYPE_PROP(Co2InjectionBaseProblem, Grid, Dune::YaspGrid<2>);
92 SET_PROP(Co2InjectionBaseProblem, FluidSystem)
96 typedef Ewoms::Co2Injection::CO2Tables CO2Tables;
99 typedef Opm::FluidSystems::BrineCO2<Scalar, CO2Tables> type;
104 SET_PROP(Co2InjectionBaseProblem, MaterialLaw)
107 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
108 enum { liquidPhaseIdx = FluidSystem::liquidPhaseIdx };
109 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
112 typedef Opm::TwoPhaseMaterialTraits<Scalar,
113 FluidSystem::liquidPhaseIdx,
114 FluidSystem::gasPhaseIdx> Traits;
118 typedef Opm::RegularizedBrooksCorey<Traits> EffMaterialLaw;
122 typedef Opm::EffToAbsLaw<EffMaterialLaw> type;
126 SET_PROP(Co2InjectionBaseProblem, HeatConductionLaw)
130 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
134 typedef Opm::Somerton<FluidSystem, Scalar> type;
138 SET_TAG_PROP(Co2InjectionBaseProblem, LinearSolverSplice, ParallelAmgLinearSolver);
141 SET_BOOL_PROP(Co2InjectionBaseProblem, NewtonWriteConvergence,
false);
148 SET_SCALAR_PROP(Co2InjectionBaseProblem, FluidSystemPressureHigh, 4e7);
149 SET_INT_PROP(Co2InjectionBaseProblem, FluidSystemNumPressure, 100);
150 SET_SCALAR_PROP(Co2InjectionBaseProblem, FluidSystemTemperatureLow, 290);
151 SET_SCALAR_PROP(Co2InjectionBaseProblem, FluidSystemTemperatureHigh, 500);
152 SET_INT_PROP(Co2InjectionBaseProblem, FluidSystemNumTemperature, 100);
156 SET_STRING_PROP(Co2InjectionBaseProblem, SimulationName,
"co2injection");
165 SET_STRING_PROP(Co2InjectionBaseProblem, GridFile,
"data/co2injection.dgf");
192 template <
class TypeTag>
195 typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
197 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
198 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
199 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
200 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
202 enum { dim = GridView::dimension };
203 enum { dimWorld = GridView::dimensionworld };
206 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
207 enum { numPhases = FluidSystem::numPhases };
208 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
209 enum { liquidPhaseIdx = FluidSystem::liquidPhaseIdx };
210 enum { CO2Idx = FluidSystem::CO2Idx };
211 enum { BrineIdx = FluidSystem::BrineIdx };
212 enum { conti0EqIdx = Indices::conti0EqIdx };
213 enum { contiCO2EqIdx = conti0EqIdx + CO2Idx };
215 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
216 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
217 typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
218 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
220 typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
221 typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
222 typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLaw) HeatConductionLaw;
223 typedef typename HeatConductionLaw::Params HeatConductionLawParams;
225 typedef Opm::MathToolbox<Evaluation> Toolbox;
226 typedef typename GridView::ctype CoordScalar;
227 typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
228 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
235 : ParentType(simulator)
243 ParentType::finishInit();
247 temperatureLow_ =
EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureLow);
248 temperatureHigh_ =
EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemTemperatureHigh);
249 nTemperature_ =
EWOMS_GET_PARAM(TypeTag,
unsigned, FluidSystemNumTemperature);
251 pressureLow_ =
EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemPressureLow);
252 pressureHigh_ =
EWOMS_GET_PARAM(TypeTag, Scalar, FluidSystemPressureHigh);
253 nPressure_ =
EWOMS_GET_PARAM(TypeTag,
unsigned, FluidSystemNumPressure);
260 FluidSystem::init(temperatureLow_,
267 fineLayerBottom_ = 22.0;
270 fineK_ = this->toDimMatrix_(1e-13);
271 coarseK_ = this->toDimMatrix_(1e-12);
275 coarsePorosity_ = 0.3;
278 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
279 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
280 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
281 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
284 fineMaterialParams_.setEntryPressure(1e4);
285 coarseMaterialParams_.setEntryPressure(5e3);
286 fineMaterialParams_.setLambda(2.0);
287 coarseMaterialParams_.setLambda(2.0);
289 fineMaterialParams_.finalize();
290 coarseMaterialParams_.finalize();
293 computeHeatCondParams_(fineHeatCondParams_, finePorosity_);
294 computeHeatCondParams_(coarseHeatCondParams_, coarsePorosity_);
302 ParentType::registerParameters();
305 "The lower temperature [K] for tabulation of the " 308 "The upper temperature [K] for tabulation of the " 311 "The number of intervals between the lower and " 312 "upper temperature");
315 "The lower pressure [Pa] for tabulation of the " 318 "The upper pressure [Pa] for tabulation of the " 321 "The number of intervals between the lower and " 325 "The temperature [K] in the reservoir");
327 "The maximum depth [m] of the reservoir");
329 "The name of the simulation used for the output " 343 std::ostringstream oss;
345 <<
"_" << Model::name();
348 oss <<
"_" << Model::discretizationName();
358 Scalar tol = this->model().newtonMethod().tolerance()*1e5;
359 this->model().checkConservativeness(tol);
362 PrimaryVariables storageL, storageG;
363 this->model().globalPhaseStorage(storageL, 0);
364 this->model().globalPhaseStorage(storageG, 1);
367 if (this->gridView().comm().rank() == 0) {
368 std::cout <<
"Storage: liquid=[" << storageL <<
"]" 369 <<
" gas=[" << storageG <<
"]\n" << std::flush;
377 template <
class Context>
378 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const 380 const auto& pos = context.pos(spaceIdx, timeIdx);
381 if (inHighTemperatureRegion_(pos))
382 return temperature_ + 100;
389 template <
class Context>
391 unsigned timeIdx)
const 393 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
394 if (isFineMaterial_(pos))
402 template <
class Context>
403 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const 405 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
406 if (isFineMaterial_(pos))
407 return finePorosity_;
408 return coarsePorosity_;
414 template <
class Context>
416 unsigned spaceIdx,
unsigned timeIdx)
const 418 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
419 if (isFineMaterial_(pos))
420 return fineMaterialParams_;
421 return coarseMaterialParams_;
429 template <
class Context>
431 unsigned spaceIdx OPM_UNUSED,
432 unsigned timeIdx OPM_UNUSED)
const 441 template <
class Context>
442 const HeatConductionLawParams &
445 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
446 if (isFineMaterial_(pos))
447 return fineHeatCondParams_;
448 return coarseHeatCondParams_;
461 template <
class Context>
462 void boundary(BoundaryRateVector& values,
const Context& context,
463 unsigned spaceIdx,
unsigned timeIdx)
const 465 const auto& pos = context.pos(spaceIdx, timeIdx);
466 if (onLeftBoundary_(pos)) {
467 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
468 initialFluidState_(fs, context, spaceIdx, timeIdx);
472 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
474 else if (onInlet_(pos)) {
475 RateVector massRate(0.0);
476 massRate[contiCO2EqIdx] = -1e-3;
478 typedef Opm::ImmiscibleFluidState<Scalar, FluidSystem> FluidState;
480 fs.setSaturation(gasPhaseIdx, 1.0);
482 context.intensiveQuantities(spaceIdx, timeIdx).fluidState().pressure(gasPhaseIdx);
483 fs.setPressure(gasPhaseIdx, Toolbox::value(pg));
484 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
486 typename FluidSystem::template ParameterCache<Scalar> paramCache;
487 paramCache.updatePhase(fs, gasPhaseIdx);
488 Scalar h = FluidSystem::template enthalpy<FluidState, Scalar>(fs, paramCache, gasPhaseIdx);
491 values.setMassRate(massRate);
492 values.setEnthalpyRate(massRate[contiCO2EqIdx] * h);
509 template <
class Context>
510 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
511 unsigned timeIdx)
const 513 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
514 initialFluidState_(fs, context, spaceIdx, timeIdx);
519 values.assignNaive(fs);
528 template <
class Context>
530 const Context& context OPM_UNUSED,
531 unsigned spaceIdx OPM_UNUSED,
532 unsigned timeIdx OPM_UNUSED)
const 533 { rate = Scalar(0.0); }
538 template <
class Context,
class Flu
idState>
539 void initialFluidState_(FluidState& fs,
540 const Context& context,
542 unsigned timeIdx)
const 544 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
549 fs.setTemperature(
temperature(context, spaceIdx, timeIdx));
554 fs.setSaturation(FluidSystem::liquidPhaseIdx, 1.0);
555 fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
560 Scalar densityL = FluidSystem::Brine::liquidDensity(temperature_, Scalar(1e5));
561 Scalar depth = maxDepth_ - pos[dim - 1];
562 Scalar pl = 1e5 - densityL * this->gravity()[dim - 1] * depth;
564 Scalar pC[numPhases];
566 MaterialLaw::capillaryPressures(pC, matParams, fs);
568 fs.setPressure(liquidPhaseIdx, pl + (pC[liquidPhaseIdx] - pC[liquidPhaseIdx]));
569 fs.setPressure(gasPhaseIdx, pl + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
574 fs.setMoleFraction(liquidPhaseIdx, CO2Idx, 0.005);
575 fs.setMoleFraction(liquidPhaseIdx, BrineIdx,
576 1.0 - fs.moleFraction(liquidPhaseIdx, CO2Idx));
578 typename FluidSystem::template ParameterCache<Scalar> paramCache;
579 typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> CFRP;
580 CFRP::solve(fs, paramCache,
586 bool onLeftBoundary_(
const GlobalPosition& pos)
const 587 {
return pos[0] < eps_; }
589 bool onRightBoundary_(
const GlobalPosition& pos)
const 590 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
592 bool onInlet_(
const GlobalPosition& pos)
const 593 {
return onRightBoundary_(pos) && (5 < pos[1]) && (pos[1] < 15); }
595 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const 596 {
return (pos[0] > 20) && (pos[0] < 30) && (pos[1] > 5) && (pos[1] < 35); }
598 void computeHeatCondParams_(HeatConductionLawParams& params, Scalar poro)
600 Scalar lambdaWater = 0.6;
601 Scalar lambdaGranite = 2.8;
603 Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
604 * std::pow(lambdaWater, poro);
605 Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
607 params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
608 params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
609 params.setVacuumLambda(lambdaDry);
612 bool isFineMaterial_(
const GlobalPosition& pos)
const 613 {
return pos[dim - 1] > fineLayerBottom_; }
617 Scalar fineLayerBottom_;
619 Scalar finePorosity_;
620 Scalar coarsePorosity_;
622 MaterialLawParams fineMaterialParams_;
623 MaterialLawParams coarseMaterialParams_;
625 HeatConductionLawParams fineHeatCondParams_;
626 HeatConductionLawParams coarseHeatCondParams_;
632 unsigned nTemperature_;
635 Scalar pressureLow_, pressureHigh_;
636 Scalar temperatureLow_, temperatureHigh_;
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
const HeatConductionLawParams & heatConductionParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:443
Scalar heatCapacitySolid(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: co2injectionproblem.hh:430
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:390
Definition: baseauxiliarymodule.hh:37
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:403
void endTimeStep()
Called by the simulator after each time integration.
Definition: co2injectionproblem.hh:355
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: co2injectionproblem.hh:241
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:415
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
#define SET_TAG_PROP(EffTypeTagName, PropTagName, ValueTypeTagName)
Define a property containing a type tag.
Definition: propertysystem.hh:436
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
Problem where is injected under a low permeable layer at a depth of 2700m.
Definition: co2injectionproblem.hh:193
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: co2injectionproblem.hh:529
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: co2injectionproblem.hh:510
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
static void registerParameters()
Definition: co2injectionproblem.hh:300
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: co2injectionproblem.hh:462
std::string name() const
The problem name.
Definition: co2injectionproblem.hh:341
Provides a linear solver backend using the parallel algebraic multi-grid (AMG) linear solver from DUN...
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
Co2InjectionProblem(Simulator &simulator)
Definition: co2injectionproblem.hh:234
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: co2injectionproblem.hh:378
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
#define SET_STRING_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant string value.
Definition: propertysystem.hh:416
A fully-implicit multi-phase flow model which assumes immiscibility of the phases.
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394