28 #ifndef EWOMS_ENERGY_MODULE_HH
29 #define EWOMS_ENERGY_MODULE_HH
34 #include <dune/common/fvector.hh>
39 namespace Properties {
53 template <
class TypeTag,
bool enableEnergy>
59 template <
class TypeTag>
63 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
64 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
65 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
66 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
67 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
68 typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
70 typedef typename FluidSystem::ParameterCache ParameterCache;
74 typedef Dune::FieldVector<Evaluation, numEq> EvalEqVector;
119 template <
class Flu
idState>
121 const FluidState &fs)
128 template <
class Flu
idState>
130 const FluidState &fluidState,
132 const Evaluation& volume)
139 const Evaluation& rate)
146 const Evaluation& rate)
154 typedef Opm::MathToolbox<Evaluation> Toolbox;
156 return Toolbox::createConstant(0.0);
163 template <
class LhsEval>
165 const IntensiveQuantities &intQuants,
173 template <
class LhsEval,
class Scv>
175 const IntensiveQuantities& intQuants,
184 template <
class LhsEval>
186 const IntensiveQuantities& intQuants)
195 template <
class Context>
197 const Context &context,
207 template <
class Context>
209 const Context &context,
220 template <
class Context>
222 const Context &context,
231 template <
class TypeTag>
235 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
236 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
238 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
239 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
240 typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
241 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
244 typedef typename FluidSystem::ParameterCache ParameterCache;
247 enum { numPhases = FluidSystem::numPhases };
248 enum { energyEqIdx = Indices::energyEqIdx };
249 enum { temperatureIdx = Indices::temperatureIdx };
251 typedef Dune::FieldVector<Evaluation, numEq> EvalEqVector;
252 typedef Opm::MathToolbox<Evaluation> Toolbox;
268 if (pvIdx == temperatureIdx)
269 return "temperature";
280 if (eqIdx == energyEqIdx)
291 if (pvIdx != temperatureIdx)
295 return std::max(1.0/1000, 1.0/model.solution(0)[globalDofIdx][temperatureIdx]);
305 if (eqIdx != energyEqIdx)
309 return 1.0 / 1.0035e3;
316 { rateVec[energyEqIdx] = rate; }
322 { rateVec[energyEqIdx] += rate; }
329 {
return -extQuants.temperatureGradNormal() * extQuants.heatConductivity(); }
335 template <
class Flu
idState>
337 const FluidState &fluidState,
int phaseIdx,
340 rateVec[energyEqIdx] =
342 * fluidState.density(phaseIdx)
343 * fluidState.enthalpy(phaseIdx);
349 template <
class Flu
idState>
351 const FluidState &fs)
353 priVars[temperatureIdx] = Toolbox::value(fs.temperature(0));
355 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
356 assert(std::abs(Toolbox::value(fs.temperature(0))
357 - Toolbox::value(fs.temperature(phaseIdx))) < 1e-30);
366 template <
class LhsEval>
368 const IntensiveQuantities &intQuants,
int phaseIdx)
370 const auto &fs = intQuants.fluidState();
371 storage[energyEqIdx] +=
372 Toolbox::template toLhs<LhsEval>(fs.density(phaseIdx))
373 * Toolbox::template toLhs<LhsEval>(fs.internalEnergy(phaseIdx))
374 * Toolbox::template toLhs<LhsEval>(fs.saturation(phaseIdx))
375 * Toolbox::template toLhs<LhsEval>(intQuants.porosity());
382 template <
class Scv,
class LhsEval>
384 const IntensiveQuantities &intQuants,
385 const Scv &scv,
int phaseIdx)
387 const auto &fs = intQuants.fractureFluidState();
388 storage[energyEqIdx] +=
389 Toolbox::template toLhs<LhsEval>(fs.density(phaseIdx))
390 * Toolbox::template toLhs<LhsEval>(fs.internalEnergy(phaseIdx))
391 * Toolbox::template toLhs<LhsEval>(fs.saturation(phaseIdx))
392 * Toolbox::template toLhs<LhsEval>(intQuants.fracturePorosity())
393 * Toolbox::template toLhs<LhsEval>(intQuants.fractureVolume())/scv.volume();
400 template <
class LhsEval>
402 const IntensiveQuantities &intQuants)
404 storage[energyEqIdx] +=
405 Toolbox::template toLhs<LhsEval>(intQuants.heatCapacitySolid())
406 * Toolbox::template toLhs<LhsEval>(intQuants.fluidState().temperature(0));
415 template <
class Context>
416 static void addAdvectiveFlux(RateVector &flux,
const Context &context,
int spaceIdx,
int timeIdx)
418 const auto &extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
421 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
422 if (!context.model().phaseIsConsidered(phaseIdx))
426 const IntensiveQuantities &up =
427 context.intensiveQuantities(extQuants.upstreamIndex(phaseIdx), timeIdx);
430 extQuants.volumeFlux(phaseIdx)
431 * up.fluidState().enthalpy(phaseIdx)
432 * up.fluidState().density(phaseIdx);
441 template <
class Context>
443 int spaceIdx,
int timeIdx)
445 const auto &scvf = context.stencil(timeIdx).interiorFace(spaceIdx);
446 const auto &extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
451 1 - extQuants.fractureWidth()/(2*scvf.area());
454 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
455 if (!context.model().phaseIsConsidered(phaseIdx))
459 const IntensiveQuantities &up =
460 context.intensiveQuantities(extQuants.upstreamIndex(phaseIdx), timeIdx);
463 extQuants.fractureVolumeFlux(phaseIdx)
464 * up.fluidState().enthalpy(phaseIdx)
465 * up.fluidState().density(phaseIdx);
475 template <
class Context>
477 int spaceIdx,
int timeIdx)
479 const auto &extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
483 - extQuants.temperatureGradNormal()
484 * extQuants.heatConductivity();
494 template <
int PVOffset,
bool enableEnergy>
500 template <
int PVOffset>
510 template <
int PVOffset>
514 enum { temperatureIdx = PVOffset };
517 enum { energyEqIdx = PVOffset };
529 template <
class TypeTag,
bool enableEnergy>
535 template <
class TypeTag>
539 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
540 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
541 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
542 typedef typename FluidSystem::ParameterCache ParameterCache;
544 typedef Opm::MathToolbox<Evaluation> Toolbox;
554 OPM_THROW(std::logic_error,
"Method heatCapacitySolid() does not make "
555 "sense for isothermal models");
565 OPM_THROW(std::logic_error,
"Method heatConductivity() does not make "
566 "sense for isothermal models");
573 template <
class Flu
idState,
class Context>
575 const Context &context,
int spaceIdx,
578 Scalar T = context.problem().temperature(context, spaceIdx, timeIdx);
579 fluidState.setTemperature(Toolbox::createConstant(T));
586 template <
class Flu
idState>
588 ParameterCache ¶mCache,
589 const ElementContext &elemCtx,
598 template <
class TypeTag>
602 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
603 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
604 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
605 typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLaw) HeatConductionLaw;
607 typedef typename FluidSystem::ParameterCache ParameterCache;
609 enum { numPhases = FluidSystem::numPhases };
610 enum { energyEqIdx = Indices::energyEqIdx };
611 enum { temperatureIdx = Indices::temperatureIdx };
613 typedef Opm::MathToolbox<Evaluation> Toolbox;
619 template <
class Flu
idState,
class Context>
621 const Context &context,
int spaceIdx,
624 const auto& priVars = context.primaryVars(spaceIdx, timeIdx);
627 val = Toolbox::createVariable(priVars[temperatureIdx], temperatureIdx);
629 val = Toolbox::createConstant(priVars[temperatureIdx]);
631 fluidState.setTemperature(val);
638 template <
class Flu
idState>
640 ParameterCache ¶mCache,
641 const ElementContext &elemCtx,
646 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
647 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
650 fs.setEnthalpy(phaseIdx,
651 FluidSystem::enthalpy(fs, paramCache, phaseIdx));
655 const auto &problem = elemCtx.problem();
656 const auto &heatCondParams = problem.heatConductionParams(elemCtx, dofIdx, timeIdx);
658 heatCapacitySolid_ = problem.heatCapacitySolid(elemCtx, dofIdx, timeIdx);
660 HeatConductionLaw::template heatConductivity<FluidState, Evaluation>(heatCondParams, fs);
662 Valgrind::CheckDefined(heatCapacitySolid_);
663 Valgrind::CheckDefined(heatConductivity_);
673 {
return heatCapacitySolid_; }
681 {
return heatConductivity_; }
684 Evaluation heatCapacitySolid_;
685 Evaluation heatConductivity_;
694 template <
class TypeTag,
bool enableEnergy>
700 template <
class TypeTag>
704 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
711 void update_(
const ElementContext &elemCtx,
int faceIdx,
int timeIdx)
714 template <
class Context,
class Flu
idState>
716 const FluidState &fs)
725 OPM_THROW(std::logic_error,
"Method temperatureGradNormal() does not "
726 "make sense for isothermal models");
735 OPM_THROW(std::logic_error,
"Method heatConductivity() does not make "
736 "sense for isothermal models");
743 template <
class TypeTag>
746 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
748 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
751 enum { dimWorld = GridView::dimensionworld };
752 typedef Dune::FieldVector<Evaluation, dimWorld> EvalDimVector;
753 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
760 void update_(
const ElementContext &elemCtx,
int faceIdx,
int timeIdx)
762 const auto& gradCalc = elemCtx.gradientCalculator();
765 EvalDimVector temperatureGrad;
766 gradCalc.calculateGradient(temperatureGrad,
769 temperatureCallback);
772 const auto &face = elemCtx.stencil(0).interiorFace(faceIdx);
774 temperatureGradNormal_ = 0;
775 for (
int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
776 temperatureGradNormal_ += (face.normal()[dimIdx]*temperatureGrad[dimIdx]);
778 const auto &extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
779 const auto &intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
780 const auto &intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
784 0.5 * (intQuantsInside.heatConductivity() + intQuantsOutside.heatConductivity());
785 Valgrind::CheckDefined(heatConductivity_);
788 template <
class Context,
class Flu
idState>
789 void updateBoundary_(
const Context &context,
int bfIdx,
int timeIdx,
const FluidState &fs)
791 const auto &stencil = context.stencil(timeIdx);
792 const auto &face = stencil.boundaryFace(bfIdx);
794 const auto &elemCtx = context.elementContext();
795 int insideScvIdx = face.interiorIndex();
796 const auto &insideScv = elemCtx.stencil(timeIdx).subControlVolume(insideScvIdx);
798 const auto &intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx);
799 const auto &fsInside = intQuantsInside.fluidState();
802 DimVector distVec = face.integrationPos();
803 distVec -= insideScv.geometry().center();
806 for (
int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
807 tmp += distVec[dimIdx] * face.normal()[dimIdx];
816 temperatureGradNormal_ =
817 (fs.temperature(0) - fsInside.temperature(0)) / dist;
820 heatConductivity_ = intQuantsInside.heatConductivity();
828 {
return temperatureGradNormal_; }
835 {
return heatConductivity_; }
838 Evaluation temperatureGradNormal_;
839 Evaluation heatConductivity_;
static void addAdvectiveFlux(RateVector &flux, const Context &context, int spaceIdx, int timeIdx)
Evaluates the advective energy fluxes for a flux integration point and adds the result in the flux ve...
Definition: energymodule.hh:416
Provides the volumetric quantities required for the energy equation.
Definition: energymodule.hh:530
static std::string eqName(int eqIdx)
Returns the name of an equation or an empty string if the specified equation index does not belong to...
Definition: energymodule.hh:278
static void addDiffusiveFlux(RateVector &flux, const Context &context, int spaceIdx, int timeIdx)
Adds the diffusive heat flux to the flux vector over the face of a sub-control volume.
Definition: energymodule.hh:476
static Evaluation heatConductionRate(const ExtensiveQuantities &extQuants)
Add the rate of the conductive heat flux to a rate vector.
Definition: energymodule.hh:152
static Scalar primaryVarWeight(const Model &model, int globalDofIdx, int pvIdx)
Returns the relative weight of a primary variable for calculating relative errors.
Definition: energymodule.hh:289
Provides the quantities required to calculate energy fluxes.
Definition: energymodule.hh:695
Scalar temperatureGradNormal() const
The temperature gradient times the face normal [K m^2 / m].
Definition: energymodule.hh:723
static void addToEnthalpyRate(RateVector &rateVec, const Evaluation &rate)
Add the rate of the enthalpy flux to a rate vector.
Definition: energymodule.hh:145
static void addFracturePhaseStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants, const Scv &scv, int phaseIdx)
Add the energy storage term for a fluid phase to an equation vector.
Definition: energymodule.hh:174
static void registerParameters()
Register all run-time parameters for the energy module.
Definition: energymodule.hh:80
static void addFracturePhaseStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants, const Scv &scv, int phaseIdx)
Add the energy storage term for a fluid phase to an equation vector.
Definition: energymodule.hh:383
const Evaluation & heatCapacitySolid() const
Returns the total heat capacity of the rock matrix in the sub-control volume.
Definition: energymodule.hh:672
Declare the properties used by the infrastructure code of the finite volume discretizations.
Provides the indices required for the energy equation.
Definition: energymodule.hh:495
Callback class for temperature.
Definition: quantitycallbacks.hh:44
static void setEnthalpyRate(RateVector &rateVec, const Evaluation &rate)
Add the rate of the enthalpy flux to a rate vector.
Definition: energymodule.hh:138
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
static void handleFractureFlux(RateVector &flux, const Context &context, int spaceIdx, int timeIdx)
Evaluates the advective energy fluxes over a fracture which should be attributed to a face of a subco...
Definition: energymodule.hh:442
static void setEnthalpyRate(RateVector &rateVec, const FluidState &fluidState, int phaseIdx, Scalar volume)
Given a fluid state, set the enthalpy rate which emerges from a volumetric rate.
Definition: energymodule.hh:336
static Evaluation heatConductionRate(const ExtensiveQuantities &extQuants)
Returns the rate of the conductive heat flux for a given flux integration point.
Definition: energymodule.hh:328
static void addToEnthalpyRate(RateVector &rateVec, const Evaluation &rate)
Add the rate of the enthalpy flux to a rate vector.
Definition: energymodule.hh:321
static std::string primaryVarName(int pvIdx)
Returns the name of a primary variable or an empty string if the specified primary variable index doe...
Definition: energymodule.hh:88
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:485
static Scalar eqWeight(const Model &model, int globalDofIdx, int eqIdx)
Returns the relative weight of a equation of the residual.
Definition: energymodule.hh:111
static void setPriVarTemperatures(PrimaryVariables &priVars, const FluidState &fs)
Given a fluid state, set the temperature in the primary variables.
Definition: energymodule.hh:350
static void setPriVarTemperatures(PrimaryVariables &priVars, const FluidState &fs)
Given a fluid state, set the temperature in the primary variables.
Definition: energymodule.hh:120
void update_(const ElementContext &elemCtx, int faceIdx, int timeIdx)
Update the quantities required to calculate energy fluxes.
Definition: energymodule.hh:760
NEW_PROP_TAG(Grid)
The type of the DUNE grid.
static Scalar primaryVarWeight(const Model &model, int globalDofIdx, int pvIdx)
Returns the relative weight of a primary variable for calculating relative errors.
Definition: energymodule.hh:103
static void addSolidHeatStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Add the energy storage term for the fracture part a fluid phase to an equation vector.
Definition: energymodule.hh:185
Provides the auxiliary methods required for consideration of the energy equation. ...
Definition: energymodule.hh:54
static void registerParameters()
Register all run-time parameters for the energy module.
Definition: energymodule.hh:258
static void addAdvectiveFlux(RateVector &flux, const Context &context, int spaceIdx, int timeIdx)
Evaluates the advective energy fluxes over a face of a subcontrol volume and adds the result in the f...
Definition: energymodule.hh:196
Definition: baseauxiliarymodule.hh:35
Scalar heatConductivity() const
The total heat conductivity at the face .
Definition: energymodule.hh:733
static void addSolidHeatStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Add the energy storage term for a fluid phase to an equation vector.
Definition: energymodule.hh:401
void update_(const ElementContext &elemCtx, int faceIdx, int timeIdx)
Update the quantities required to calculate energy fluxes.
Definition: energymodule.hh:711
static void addPhaseStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants, int phaseIdx)
Add the energy storage term for a fluid phase to an equation vector.
Definition: energymodule.hh:367
static void addPhaseStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants, int phaseIdx)
Add the energy storage term for a fluid phase to an equation vector.
Definition: energymodule.hh:164
static void updateTemperatures_(FluidState &fluidState, const Context &context, int spaceIdx, int timeIdx)
Update the temperatures of the fluids of a fluid state.
Definition: energymodule.hh:620
This method contains all callback classes for quantities that are required by some extensive quantiti...
const Evaluation & temperatureGradNormal() const
The temperature gradient times the face normal [K m^2 / m].
Definition: energymodule.hh:827
void updateBoundary_(const Context &context, int bfIdx, int timeIdx, const FluidState &fs)
Definition: energymodule.hh:715
static void updateTemperatures_(FluidState &fluidState, const Context &context, int spaceIdx, int timeIdx)
Update the temperatures of the fluids of a fluid state.
Definition: energymodule.hh:574
Evaluation heatConductivity() const
Returns the total conductivity capacity of the rock matrix in the sub-control volume.
Definition: energymodule.hh:563
void update_(FluidState &fs, ParameterCache ¶mCache, const ElementContext &elemCtx, int dofIdx, int timeIdx)
Update the quantities required to calculate energy fluxes.
Definition: energymodule.hh:587
static void setEnthalpyRate(RateVector &rateVec, const Evaluation &rate)
Add the rate of the enthalpy flux to a rate vector.
Definition: energymodule.hh:315
static std::string primaryVarName(int pvIdx)
Returns the name of a primary variable or an empty string if the specified primary variable index doe...
Definition: energymodule.hh:266
const Evaluation & heatConductivity() const
The total heat conductivity at the face .
Definition: energymodule.hh:834
void update_(FluidState &fs, ParameterCache ¶mCache, const ElementContext &elemCtx, int dofIdx, int timeIdx)
Update the quantities required to calculate energy fluxes.
Definition: energymodule.hh:639
static void setEnthalpyRate(RateVector &rateVec, const FluidState &fluidState, int phaseIdx, const Evaluation &volume)
Given a fluid state, set the enthalpy rate which emerges from a volumetric rate.
Definition: energymodule.hh:129
Evaluation heatCapacitySolid() const
Returns the total heat capacity of the rock matrix in the sub-control volume.
Definition: energymodule.hh:552
static std::string eqName(int eqIdx)
Returns the name of an equation or an empty string if the specified equation index does not belong to...
Definition: energymodule.hh:96
static void handleFractureFlux(RateVector &flux, const Context &context, int spaceIdx, int timeIdx)
Evaluates the advective energy fluxes over a fracture which should be attributed to a face of a subco...
Definition: energymodule.hh:208
const Evaluation & heatConductivity() const
Returns the total conductivity capacity of the rock matrix in the sub-control volume.
Definition: energymodule.hh:680
static void addDiffusiveFlux(RateVector &flux, const Context &context, int spaceIdx, int timeIdx)
Adds the diffusive heat flux to the flux vector over the face of a sub-control volume.
Definition: energymodule.hh:221
void updateBoundary_(const Context &context, int bfIdx, int timeIdx, const FluidState &fs)
Definition: energymodule.hh:789
static Scalar eqWeight(const Model &model, int globalDofIdx, int eqIdx)
Returns the relative weight of a equation.
Definition: energymodule.hh:301