24#ifndef DUMUX_GRAVITYPART_HH
25#define DUMUX_GRAVITYPART_HH
47template<
class TypeTag>
65 dim = GridView::dimension, dimWorld = GridView::dimensionworld
69 wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, numPhases = getPropValue<TypeTag, Properties::NumPhases>()
72 using Intersection =
typename GridView::Intersection;
73 using DimVector = Dune::FieldVector<Scalar, dim>;
74 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
75 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
87 void getFlux(DimVector& flux,
const Intersection& intersection,
const Scalar satI,
const Scalar satJ)
const
89 auto element = intersection.inside();
91 int globalIdxI = problem_.variables().index(element);
92 CellData& cellDataI = problem_.variables().cellData(globalIdxI);
94 int indexInInside = intersection.indexInInside();
100 Scalar lambdaNwJ = 0;
102 const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(element.geometry().center());
105 lambdaWI=cellDataI.mobility(wPhaseIdx);
106 lambdaNwI=cellDataI.mobility(nPhaseIdx);
110 lambdaWI = fluidMatrixInteraction.krw(satI);
111 lambdaWI /= viscosity_[wPhaseIdx];
112 lambdaNwI = fluidMatrixInteraction.krn(satI);
113 lambdaNwI /= viscosity_[nPhaseIdx];
116 Scalar potentialDiffW = cellDataI.fluxData().upwindPotential(wPhaseIdx, indexInInside);
117 Scalar potentialDiffNw = cellDataI.fluxData().upwindPotential(nPhaseIdx, indexInInside);
119 DimMatrix meanPermeability(0);
120 GlobalPosition distVec(0);
124 if (intersection.neighbor())
127 auto neighbor = intersection.outside();
129 int globalIdxJ = problem_.variables().index(neighbor);
130 CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
132 distVec = neighbor.geometry().center() - element.geometry().center();
135 problem_.spatialParams().meanK(meanPermeability,
136 problem_.spatialParams().intrinsicPermeability(element),
137 problem_.spatialParams().intrinsicPermeability(neighbor));
142 lambdaWJ=cellDataJ.mobility(wPhaseIdx);
143 lambdaNwJ=cellDataJ.mobility(nPhaseIdx);
147 const auto fluidMatrixInteractionNeighbor = problem_.spatialParams().fluidMatrixInteractionAtPos(neighbor.geometry().center());
148 lambdaWJ = fluidMatrixInteractionNeighbor.krw(satJ);
149 lambdaWJ /= viscosity_[wPhaseIdx];
150 lambdaNwJ = fluidMatrixInteractionNeighbor.krn(satJ);
151 lambdaNwJ /= viscosity_[nPhaseIdx];
154 lambdaW = (potentialDiffW >= 0) ? lambdaWI : lambdaWJ;
155 lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
156 lambdaNw = (potentialDiffNw >= 0) ? lambdaNwI : lambdaNwJ;
157 lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNw;
162 problem_.spatialParams().meanK(meanPermeability,
163 problem_.spatialParams().intrinsicPermeability(element));
165 distVec = intersection.geometry().center() - element.geometry().center();
168 lambdaWJ = fluidMatrixInteraction.krw(satJ);
169 lambdaWJ /= viscosity_[wPhaseIdx];
170 lambdaNwJ = fluidMatrixInteraction.krn(satJ);
171 lambdaNwJ /= viscosity_[nPhaseIdx];
174 lambdaW = (potentialDiffW > 0) ? lambdaWI : lambdaWJ;
175 lambdaNw = (potentialDiffNw > 0) ? lambdaNwI : lambdaNwJ;
179 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
180 Scalar dist = distVec.two_norm();
183 Scalar areaScaling = (unitOuterNormal * distVec);
190 Scalar scalarGravity = problem_.gravity() * distVec;
192 flux = unitOuterNormal;
195 flux *= lambdaW*lambdaNw/(lambdaW+lambdaNw) * scalarPerm * (density_[wPhaseIdx] - density_[nPhaseIdx]) * scalarGravity * areaScaling;
210 auto element = *problem_.gridView().template begin<0> ();
211 FluidState fluidState;
212 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
213 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
214 fluidState.setTemperature(problem_.temperature(element));
215 fluidState.setSaturation(wPhaseIdx, 1.);
216 fluidState.setSaturation(nPhaseIdx, 0.);
225 const bool preComput_;
227 Scalar density_[numPhases];
228 Scalar viscosity_[numPhases];
Base class for defining a convective part of the saturation transport equation.
constexpr auto getPropValue()
get the value data member of a property
Definition: propertysystem.hh:154
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
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:208
GravityPart(Problem &problem)
Constructs a GravityPart object.
Definition: gravitypart.hh:203
void getFlux(DimVector &flux, const Intersection &intersection, const Scalar satI, const Scalar satJ) const
Returns convective term for current element face.
Definition: gravitypart.hh:87
Base file for properties related to sequential models.