richardsintensivequantities.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef EWOMS_RICHARDS_INTENSIVE_QUANTITIES_HH
29#define EWOMS_RICHARDS_INTENSIVE_QUANTITIES_HH
30
31#include "richardsproperties.hh"
32
33#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
34
35#include <dune/common/fvector.hh>
36#include <dune/common/fmatrix.hh>
37
38namespace Opm {
39
46template <class TypeTag>
48 : public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
49 , public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
50{
59
61 enum { pressureWIdx = Indices::pressureWIdx };
62 enum { numPhases = FluidSystem::numPhases };
63 enum { liquidPhaseIdx = getPropValue<TypeTag, Properties::LiquidPhaseIndex>() };
64 enum { gasPhaseIdx = getPropValue<TypeTag, Properties::GasPhaseIndex>() };
65 enum { dimWorld = GridView::dimensionworld };
66
67 using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities;
68 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
69 using ScalarPhaseVector = Dune::FieldVector<Scalar, numPhases>;
70 using PhaseVector = Dune::FieldVector<Evaluation, numPhases>;
71 using Toolbox = Opm::MathToolbox<Evaluation>;
72
73public:
75 using FluidState = Opm::ImmiscibleFluidState<Evaluation, FluidSystem>;
76
78 {}
79
81
83
87 void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
88 {
89 ParentType::update(elemCtx, dofIdx, timeIdx);
90
91 const auto& T = elemCtx.problem().temperature(elemCtx, dofIdx, timeIdx);
92 fluidState_.setTemperature(T);
93
94 // material law parameters
95 const auto& problem = elemCtx.problem();
96 const typename MaterialLaw::Params& materialParams =
97 problem.materialLawParams(elemCtx, dofIdx, timeIdx);
98 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
99
101 // calculate the pressures
103
104 // first, we have to find the minimum capillary pressure (i.e. Sw = 0)
105 fluidState_.setSaturation(liquidPhaseIdx, 1.0);
106 fluidState_.setSaturation(gasPhaseIdx, 0.0);
107 ScalarPhaseVector pC;
108 MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
109
110 // non-wetting pressure can be larger than the
111 // reference pressure if the medium is fully
112 // saturated by the wetting phase
113 const Evaluation& pW = priVars.makeEvaluation(pressureWIdx, timeIdx);
114 Evaluation pN =
115 Toolbox::max(elemCtx.problem().referencePressure(elemCtx, dofIdx, /*timeIdx=*/0),
116 pW + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
117
119 // calculate the saturations
121 fluidState_.setPressure(liquidPhaseIdx, pW);
122 fluidState_.setPressure(gasPhaseIdx, pN);
123
124 PhaseVector sat;
125 MaterialLaw::saturations(sat, materialParams, fluidState_);
126 fluidState_.setSaturation(liquidPhaseIdx, sat[liquidPhaseIdx]);
127 fluidState_.setSaturation(gasPhaseIdx, sat[gasPhaseIdx]);
128
129 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
130 paramCache.updateAll(fluidState_);
131
132 // compute and set the wetting phase viscosity
133 const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, liquidPhaseIdx);
134 fluidState_.setViscosity(liquidPhaseIdx, mu);
135 fluidState_.setViscosity(gasPhaseIdx, 1e-20);
136
137 // compute and set the wetting phase density
138 const Evaluation& rho = FluidSystem::density(fluidState_, paramCache, liquidPhaseIdx);
139 fluidState_.setDensity(liquidPhaseIdx, rho);
140 fluidState_.setDensity(gasPhaseIdx, 1e-20);
141
142 // relperms
143 MaterialLaw::relativePermeabilities(relativePermeability_, materialParams, fluidState_);
144
145 // mobilities
146 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
147 mobility_[phaseIdx] = relativePermeability_[phaseIdx]/fluidState_.viscosity(phaseIdx);
148
149 // porosity
150 porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
151
152 // intrinsic permeability
153 intrinsicPerm_ = problem.intrinsicPermeability(elemCtx, dofIdx, timeIdx);
154
155 // update the quantities specific for the velocity model
156 FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
157 }
158
162 const FluidState& fluidState() const
163 { return fluidState_; }
164
168 const Evaluation& porosity() const
169 { return porosity_; }
170
174 const DimMatrix& intrinsicPermeability() const
175 { return intrinsicPerm_; }
176
180 const Evaluation& relativePermeability(unsigned phaseIdx) const
181 { return relativePermeability_[phaseIdx]; }
182
186 const Evaluation& mobility(unsigned phaseIdx) const
187 { return mobility_[phaseIdx]; }
188
189private:
190 FluidState fluidState_;
191 DimMatrix intrinsicPerm_;
192 std::array<Evaluation,numPhases> relativePermeability_;
193 std::array<Evaluation,numPhases> mobility_;
194 Evaluation porosity_;
195};
196
197} // namespace Opm
198
199#endif
Intensive quantities required by the Richards model.
Definition: richardsintensivequantities.hh:50
const Evaluation & mobility(unsigned phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: richardsintensivequantities.hh:186
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: richardsintensivequantities.hh:162
const Evaluation & porosity() const
Returns the average porosity within the control volume.
Definition: richardsintensivequantities.hh:168
RichardsIntensiveQuantities()
Definition: richardsintensivequantities.hh:77
const DimMatrix & intrinsicPermeability() const
Returns the intrinsic permeability tensor a degree of freedom.
Definition: richardsintensivequantities.hh:174
void update(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Definition: richardsintensivequantities.hh:87
Opm::ImmiscibleFluidState< Evaluation, FluidSystem > FluidState
The type returned by the fluidState() method.
Definition: richardsintensivequantities.hh:75
RichardsIntensiveQuantities(const RichardsIntensiveQuantities &other)=default
RichardsIntensiveQuantities & operator=(const RichardsIntensiveQuantities &other)=default
const Evaluation & relativePermeability(unsigned phaseIdx) const
Returns the relative permeability of a given phase within the control volume.
Definition: richardsintensivequantities.hh:180
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
Contains the property declarations for the Richards model.