26 #ifndef EWOMS_STOKES_BOUNDARY_RATE_VECTOR_HH
27 #define EWOMS_STOKES_BOUNDARY_RATE_VECTOR_HH
29 #include <opm/material/localad/Math.hpp>
30 #include <opm/material/common/Valgrind.hpp>
31 #include <opm/material/constraintsolvers/NcpFlash.hpp>
33 #include <dune/common/fvector.hh>
45 template <
class TypeTag>
48 typedef typename GET_PROP_TYPE(TypeTag, RateVector) ParentType;
49 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
50 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
51 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
52 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
54 enum { numComponents = FluidSystem::numComponents };
56 enum { dimWorld = GridView::dimensionworld };
58 enum { conti0EqIdx = Indices::conti0EqIdx };
59 enum { momentum0EqIdx = Indices::momentum0EqIdx };
63 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
93 template <
class Context,
class Flu
idState>
94 void setFreeFlow(
const Context &context,
int bfIdx,
int timeIdx,
95 const DimVector &velocity,
const FluidState &fluidState)
97 const auto &stencil = context.stencil(timeIdx);
98 const auto &scvf = stencil.boundaryFace(bfIdx);
100 int insideScvIdx = context.interiorScvIndex(bfIdx, timeIdx);
102 const auto &insideIntQuants = context.intensiveQuantities(bfIdx, timeIdx);
105 const auto &normal = scvf.normal();
108 DimVector distVec = stencil.subControlVolume(insideScvIdx).geometry().center();
109 const auto &scvPos = context.element().geometry().corner(insideScvIdx);
110 distVec.axpy(-1, scvPos);
113 for (
int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
114 dist += distVec[dimIdx] * normal[dimIdx];
115 dist = std::abs(dist);
117 DimVector gradv[dimWorld];
118 for (
int axisIdx = 0; axisIdx < dimWorld; ++axisIdx) {
121 gradv[axisIdx] = normal;
122 gradv[axisIdx] *= (velocity[axisIdx]
123 - insideIntQuants.velocity()[axisIdx]) / dist;
124 Valgrind::CheckDefined(gradv[axisIdx]);
128 Scalar volumeFlux = 0;
129 for (
int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
130 volumeFlux += velocity[dimIdx]*normal[dimIdx];
132 typename FluidSystem::ParameterCache paramCache;
133 paramCache.updatePhase(fluidState, phaseIdx);
134 Scalar density = FluidSystem::density(fluidState, paramCache, phaseIdx);
135 Scalar molarDensity = density / fluidState.averageMolarMass(phaseIdx);
136 for (
int compIdx = 0; compIdx < numComponents; ++compIdx) {
137 (*this)[conti0EqIdx + compIdx] =
140 * fluidState.moleFraction(phaseIdx, compIdx);
144 for (
int axisIdx = 0; axisIdx < dimWorld; ++axisIdx) {
147 for (
int j = 0; j < dimWorld; ++j) {
148 tmp[j] = gradv[axisIdx][j] + gradv[j][axisIdx];
153 for (
int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
154 tmp2 += tmp[dimIdx]*normal[dimIdx];
156 (*this)[momentum0EqIdx + axisIdx] = -insideIntQuants.fluidState().viscosity(phaseIdx) * tmp2;
159 EnergyModule::setEnthalpyRate(*
this, fluidState, phaseIdx, volumeFlux);
174 template <
class Context,
class Flu
idState>
175 void setInFlow(
const Context &context,
int bfIdx,
int timeIdx,
176 const DimVector &velocity,
const FluidState &fluidState)
178 const auto &intQuants = context.intensiveQuantities(bfIdx, timeIdx);
180 setFreeFlow(context, bfIdx, timeIdx, velocity, fluidState);
183 for (
int compIdx = 0; compIdx < numComponents; ++compIdx)
184 (*
this)[conti0EqIdx + compIdx] = std::min(0.0, (*
this)[conti0EqIdx + compIdx]);
187 for (
int axisIdx = 0; axisIdx < dimWorld; ++axisIdx)
188 (*
this)[momentum0EqIdx + axisIdx] = std::min(0.0, (*
this)[momentum0EqIdx + axisIdx]);
196 template <
class Context>
197 void setOutFlow(
const Context &context,
int spaceIdx,
int timeIdx)
199 const auto &intQuants = context.intensiveQuantities(spaceIdx, timeIdx);
201 DimVector velocity = intQuants.velocity();
202 const auto &fluidState = intQuants.fluidState();
204 setFreeFlow(context, spaceIdx, timeIdx, velocity, fluidState);
207 for (
int compIdx = 0; compIdx < numComponents; ++compIdx)
208 (*
this)[conti0EqIdx + compIdx] = std::max(0.0, (*
this)[conti0EqIdx + compIdx]);
211 for (
int axisIdx = 0; axisIdx < dimWorld; ++axisIdx)
212 (*
this)[momentum0EqIdx + axisIdx] = std::max(0.0, (*
this)[momentum0EqIdx + axisIdx]);
220 template <
class Context>
221 void setNoFlow(
const Context &context,
int spaceIdx,
int timeIdx)
223 static DimVector v0(0.0);
225 const auto &intQuants = context.intensiveQuantities(spaceIdx, timeIdx);
226 const auto &fluidState = intQuants.fluidState();
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:485
void setInFlow(const Context &context, int bfIdx, int timeIdx, const DimVector &velocity, const FluidState &fluidState)
Set a in-flow boundary in the (Navier-)Stoke model.
Definition: stokesboundaryratevector.hh:175
void setNoFlow(const Context &context, int spaceIdx, int timeIdx)
Set a no-flow boundary in the (Navier-)Stoke model.
Definition: stokesboundaryratevector.hh:221
void setFreeFlow(const Context &context, int bfIdx, int timeIdx, const DimVector &velocity, const FluidState &fluidState)
Definition: stokesboundaryratevector.hh:94
Provides the auxiliary methods required for consideration of the energy equation. ...
Definition: energymodule.hh:54
void setOutFlow(const Context &context, int spaceIdx, int timeIdx)
Set a out-flow boundary in the (Navier-)Stoke model.
Definition: stokesboundaryratevector.hh:197
Definition: baseauxiliarymodule.hh:35
StokesBoundaryRateVector(const StokesBoundaryRateVector &value)
Definition: stokesboundaryratevector.hh:80
Contains the intensive quantities of the Stokes model.
StokesBoundaryRateVector(Scalar value)
Definition: stokesboundaryratevector.hh:73
StokesBoundaryRateVector()
Definition: stokesboundaryratevector.hh:66
Implements a boundary vector for the fully implicit (Navier-)Stokes model.
Definition: stokesboundaryratevector.hh:46