3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
gravitypart.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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program 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 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_GRAVITYPART_HH
25#define DUMUX_GRAVITYPART_HH
26
28#include "properties.hh"
29
30namespace Dumux {
47template<class TypeTag>
48class GravityPart: public ConvectivePart<TypeTag>
49{
50private:
51 using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
52 using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
53 using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
54 using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
55
56 using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
57 using MaterialLaw = typename SpatialParams::MaterialLaw;
58
59 using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
60 using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
61
62 using CellData = typename GET_PROP_TYPE(TypeTag, CellData);
63
64 enum
65 {
66 dim = GridView::dimension, dimWorld = GridView::dimensionworld
67 };
68 enum
69 {
70 wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
71 };
72
73 using Intersection = typename GridView::Intersection;
74 using DimVector = Dune::FieldVector<Scalar, dim>;
75 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
76 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
77
78public:
79
88 void getFlux(DimVector& flux, const Intersection& intersection, const Scalar satI, const Scalar satJ) const
89 {
90 auto element = intersection.inside();
91
92 int globalIdxI = problem_.variables().index(element);
93 CellData& cellDataI = problem_.variables().cellData(globalIdxI);
94
95 int indexInInside = intersection.indexInInside();
96
97 //get lambda_bar = lambda_n*f_w
98 Scalar lambdaWI = 0;
99 Scalar lambdaNwI = 0;
100 Scalar lambdaWJ = 0;
101 Scalar lambdaNwJ = 0;
102
103 if (preComput_)
104 {
105 lambdaWI=cellDataI.mobility(wPhaseIdx);
106 lambdaNwI=cellDataI.mobility(nPhaseIdx);
107 }
108 else
109 {
110 lambdaWI = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satI);
111 lambdaWI /= viscosity_[wPhaseIdx];
112 lambdaNwI = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satI);
113 lambdaNwI /= viscosity_[nPhaseIdx];
114 }
115
116 Scalar potentialDiffW = cellDataI.fluxData().upwindPotential(wPhaseIdx, indexInInside);
117 Scalar potentialDiffNw = cellDataI.fluxData().upwindPotential(nPhaseIdx, indexInInside);
118
119 DimMatrix meanPermeability(0);
120 GlobalPosition distVec(0);
121 Scalar lambdaW = 0;
122 Scalar lambdaNw = 0;
123
124 if (intersection.neighbor())
125 {
126 // access neighbor
127 auto neighbor = intersection.outside();
128
129 int globalIdxJ = problem_.variables().index(neighbor);
130 CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
131
132 distVec = neighbor.geometry().center() - element.geometry().center();
133
134 // get permeability
135 problem_.spatialParams().meanK(meanPermeability,
136 problem_.spatialParams().intrinsicPermeability(element),
137 problem_.spatialParams().intrinsicPermeability(neighbor));
138
139 //get lambda_bar = lambda_n*f_w
140 if (preComput_)
141 {
142 lambdaWJ=cellDataJ.mobility(wPhaseIdx);
143 lambdaNwJ=cellDataJ.mobility(nPhaseIdx);
144 }
145 else
146 {
147 lambdaWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(neighbor), satJ);
148 lambdaWJ /= viscosity_[wPhaseIdx];
149 lambdaNwJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(neighbor), satJ);
150 lambdaNwJ /= viscosity_[nPhaseIdx];
151 }
152
153 lambdaW = (potentialDiffW >= 0) ? lambdaWI : lambdaWJ;
154 lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
155 lambdaNw = (potentialDiffNw >= 0) ? lambdaNwI : lambdaNwJ;
156 lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNw;
157 }
158 else
159 {
160 // get permeability
161 problem_.spatialParams().meanK(meanPermeability,
162 problem_.spatialParams().intrinsicPermeability(element));
163
164 distVec = intersection.geometry().center() - element.geometry().center();
165
166 //calculate lambda_n*f_w at the boundary
167 lambdaWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satJ);
168 lambdaWJ /= viscosity_[wPhaseIdx];
169 lambdaNwJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satJ);
170 lambdaNwJ /= viscosity_[nPhaseIdx];
171
172 //If potential is zero always take value from the boundary!
173 lambdaW = (potentialDiffW > 0) ? lambdaWI : lambdaWJ;
174 lambdaNw = (potentialDiffNw > 0) ? lambdaNwI : lambdaNwJ;
175 }
176
177 // set result to K*grad(pc)
178 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
179 Scalar dist = distVec.two_norm();
180 //calculate unit distVec
181 distVec /= dist;
182 Scalar areaScaling = (unitOuterNormal * distVec);
183
184 Dune::FieldVector<Scalar, dim> permeability(0);
185 meanPermeability.mv(unitOuterNormal, permeability);
186
187 Scalar scalarPerm = permeability.two_norm();
188
189 Scalar scalarGravity = problem_.gravity() * distVec;
190
191 flux = unitOuterNormal;
192
193 // set result to f_w*lambda_n*K*grad(pc)
194 flux *= lambdaW*lambdaNw/(lambdaW+lambdaNw) * scalarPerm * (density_[wPhaseIdx] - density_[nPhaseIdx]) * scalarGravity * areaScaling;
195 }
196
202 GravityPart (Problem& problem)
203 : ConvectivePart<TypeTag>(problem), problem_(problem), preComput_(GET_PROP_VALUE(TypeTag, PrecomputedConstRels))
204 {}
205
208 {
209 auto element = *problem_.gridView().template begin<0> ();
210 FluidState fluidState;
211 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
212 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
213 fluidState.setTemperature(problem_.temperature(element));
214 fluidState.setSaturation(wPhaseIdx, 1.);
215 fluidState.setSaturation(nPhaseIdx, 0.);
216 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
217 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
218 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
219 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
220 }
221
222private:
223 Problem& problem_; //problem data
224 const bool preComput_;//if preCompute = true the mobilities are taken from the variable object,
225 //if preCompute = false new mobilities will be taken (for implicit Scheme)
226 Scalar density_[numPhases];
227 Scalar viscosity_[numPhases];
228};
229} // end namespace Dumux
230
231#endif
#define GET_PROP_VALUE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:282
#define GET_PROP_TYPE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:283
Base class for defining a convective part of the saturation transport equation.
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Property tag PrecomputedConstRels
Definition: porousmediumflow/2p/sequential/transport/cellcentered/properties.hh:46
Property tag NumPhases
Number of phases in the system.
Definition: porousmediumflow/sequential/properties.hh:69
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Base class for defining a convective part of the saturation transport equation.
Definition: convectivepart.hh:39
Class for defining the gravity term of a two-phase flow saturation equation.
Definition: gravitypart.hh:49
void initialize()
For initialization.
Definition: gravitypart.hh:207
GravityPart(Problem &problem)
Constructs a GravityPart object.
Definition: gravitypart.hh:202
void getFlux(DimVector &flux, const Intersection &intersection, const Scalar satI, const Scalar satJ) const
Returns convective term for current element face.
Definition: gravitypart.hh:88
Base file for properties related to sequential models.