Commit cb836d2b authored by Klaus Mosthaf's avatar Klaus Mosthaf
Browse files

Enabled outflow conditions (beta version) for the 2p2cni model.

boxspatialparameters: additional empty matrixHeatFlux() method, which uses the face as argument and can handle both, boundary faces and sub-control-volume faces; this is used for the calculation of outflow fluxes only so far and has to be specified in the spatial parameters, if required. 
2p2cnifluxvariables: additional call to matrixHeatFlux(), if we are on a boundary face; changed temperature gradient to a private object with a return function to make it accessible from outside.
boxlocalresidual: removed some superfluous asImp()s.
2p2cnilocalresidual: added computeOutflowValues() method for the evaluation of outflow conditions.
regularizedvangenuchten: corrected typo.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7368 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 93ed96bf
......@@ -99,6 +99,9 @@ public:
Scalar normalMatrixHeatFlux() const
{ return normalMatrixHeatFlux_; }
Vector temperatureGradient() const
{ return temperatureGrad_; }
protected:
template<class FaceType>
void calculateValues_(const Problem &problem,
......@@ -112,29 +115,39 @@ protected:
// calculate temperature gradient using finite element
// gradients
Vector temperatureGrad(0);
temperatureGrad_ = 0;
Vector tmp(0.0);
for (int vertIdx = 0; vertIdx < this->fvGeom_.numVertices; vertIdx++)
{
tmp = this->fvGeom_.subContVolFace[scvfIdx_].grad[vertIdx];
tmp = face.grad[vertIdx];
tmp *= elemVolVars[vertIdx].temperature();
temperatureGrad += tmp;
temperatureGrad_ += tmp;
}
// The spatial parameters calculates the actual heat flux vector
problem.spatialParameters().matrixHeatFlux(tmp,
*this,
elemVolVars,
temperatureGrad,
element,
this->fvGeom_,
scvfIdx_);
if (!onBoundary)
problem.spatialParameters().matrixHeatFlux(tmp,
*this,
elemVolVars,
temperatureGrad_,
element,
this->fvGeom_,
scvfIdx_);
else
problem.spatialParameters().matrixHeatFlux(tmp,
*this,
elemVolVars,
face,
element,
this->fvGeom_);
// project the heat flux vector on the face's normal vector
normalMatrixHeatFlux_ = tmp*this->fvGeom_.subContVolFace[scvfIdx_].normal;
normalMatrixHeatFlux_ = tmp*face.normal;
}
private:
Scalar normalMatrixHeatFlux_;
Vector temperatureGrad_;
int scvfIdx_;
};
......
......@@ -52,10 +52,12 @@ class TwoPTwoCNILocalResidual : public TwoPTwoCLocalResidual<TypeTag>
typedef TwoPTwoCLocalResidual<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryVariables) BoundaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, TwoPTwoCIndices) Indices;
enum {
......@@ -174,6 +176,45 @@ public:
fluxData.normalMatrixHeatFlux();
}
/*!
* \brief Compute the fluxes at outflow boundaries. This does essentially the same
* as computeFluxes, but the fluxes are evaluated at the integration point
* of the boundary face. Therefore, still some variables are evaluated
* at the vertex (usually the ones which are upwinded).
*
* \param flux A temporary vector, where the outflow boundary fluxes are written into
* \param boundaryVars The boundary variables object
* \param scvIdx The index of the SCV containing the outflow boundary face
* \param boundaryFaceIdx The index of the boundary face
*/
void computeOutflowValues(PrimaryVariables &flux,
const BoundaryVariables &boundaryVars,
const int scvIdx,
const int boundaryFaceIdx)
{
ParentType::computeOutflowValues(flux, boundaryVars, scvIdx, boundaryFaceIdx);
const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
const VolumeVariables &vertVars = this->curVolVars_()[scvIdx];
if (bcTypes.isOutflow(energyEqIdx))
{
// advective heat flux in all phases
flux[energyEqIdx] = 0;
for (int phase = 0; phase < numPhases; ++phase) {
flux[energyEqIdx] +=
boundaryVars.KmvpNormal(phase) *
vertVars.density(phase) *
vertVars.mobility(phase) *
vertVars.enthalpy(phase);
}
// add conductive heat flux in all phases
flux[energyEqIdx] +=
boundaryVars.normalMatrixHeatFlux();
}
}
private:
Scalar massUpwindWeight_;
......
......@@ -363,10 +363,10 @@ protected:
assert(0 <= pvIdx && pvIdx < numEq);
Valgrind::CheckDefined(tmp[pvIdx]);
this->residual_[i][eqIdx] =
curPrimaryVar_(i, pvIdx) - tmp[pvIdx];
residual_[i][eqIdx] =
curPrimaryVar_(i, pvIdx) - tmp[pvIdx];
this->storageTerm_[i][eqIdx] = 0.0;
storageTerm_[i][eqIdx] = 0.0;
};
};
}
......@@ -406,15 +406,15 @@ protected:
// add the residual of all vertices of the boundary
// segment
evalNeumannSegment_(isIt,
elemVertIdx,
boundaryFaceIdx);
asImp_().evalNeumannSegment_(isIt,
elemVertIdx,
boundaryFaceIdx);
// evaluate the outflow conditions at the boundary face
// ATTENTION: This is so far a beta version that is only for the 2p2c and 2p2cni model
// available and not thoroughly tested.
this->asImp_().evalOutflowSegment(isIt,
elemVertIdx,
boundaryFaceIdx);
asImp_().evalOutflowSegment(isIt,
elemVertIdx,
boundaryFaceIdx);
}
}
}
......@@ -471,7 +471,7 @@ protected:
PrimaryVariables flux;
Valgrind::SetUndefined(flux);
this->asImp_().computeFlux(flux, k);
asImp_().computeFlux(flux, k);
Valgrind::CheckDefined(flux);
Scalar extrusionFactor =
......@@ -512,7 +512,7 @@ protected:
// all sub control volumes
for (int i=0; i < fvElemGeom_().numVertices; i++) {
Valgrind::SetUndefined(storageTerm_[i]);
this->asImp_().computeStorage(storageTerm_[i], i, /*isOldSol=*/false);
asImp_().computeStorage(storageTerm_[i], i, /*isOldSol=*/false);
storageTerm_[i] *=
fvElemGeom_().subContVol[i].volume
* curVolVars_(i).extrusionFactor();
......@@ -543,8 +543,8 @@ protected:
// doing the time discretization...
Valgrind::SetUndefined(storageTerm_[i]);
Valgrind::SetUndefined(tmp);
this->asImp_().computeStorage(storageTerm_[i], i, false);
this->asImp_().computeStorage(tmp, i, true);
asImp_().computeStorage(storageTerm_[i], i, false);
asImp_().computeStorage(tmp, i, true);
Valgrind::CheckDefined(storageTerm_[i]);
Valgrind::CheckDefined(tmp);
......@@ -557,7 +557,7 @@ protected:
// subtract the source term from the local rate
Valgrind::SetUndefined(tmp);
this->asImp_().computeSource(tmp, i);
asImp_().computeSource(tmp, i);
Valgrind::CheckDefined(tmp);
tmp *= fvElemGeom_().subContVol[i].volume * extrusionFactor;
residual_[i] -= tmp;
......
......@@ -116,7 +116,7 @@ public:
// use spline between threshold Swe and 1.0
Scalar mTh = VanGenuchten::dpC_dSw(params, SwThHigh);
Spline<Scalar> sp(SwThHigh, 1.0, // x0, x1
yTh, 0, // m0, m1
yTh, 0, // y0, y1
mTh, m1); // m0, m1
return sp.eval(Swe);
}
......
......@@ -99,6 +99,34 @@ public:
"a materialLawParamsAtPos() method.");
}
/*!
* \brief Calculate the heat flux \f$[W/m^2]\f$ through the
* rock matrix based on the temperature gradient \f$[K / m]\f$
*
* This is only required for non-isothermal models that use outflow
* boundary conditions.
*
* \param heatFlux The resulting heat flux vector
* \param fluxDat The flux variables
* \param vDat The volume variables
* \param face The boundary or SCV face
* \param element The current finite element
* \param fvElemGeom The finite volume geometry of the current element
* \tparam FaceType The type of the face (boundary face / SCV face)
*/
template <class FaceType>
void matrixHeatFlux(Vector &heatFlux,
const FluxVariables &fluxDat,
const ElementVolumeVariables &vDat,
const FaceType &face,
const Element &element,
const FVElementGeometry &fvElemGeom) const
{
DUNE_THROW(Dune::InvalidStateException,
"The spatial parameters do not provide "
"a matrixHeatFlux() method.");
}
private:
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment