diff --git a/dumux/boxmodels/2p2cni/2p2cnifluxvariables.hh b/dumux/boxmodels/2p2cni/2p2cnifluxvariables.hh
index 52d09ba6633970b6ed73de0553965a5c8c8c3463..00e62b05dca7db60e28fbbba425ec2835e18a7ce 100644
--- a/dumux/boxmodels/2p2cni/2p2cnifluxvariables.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnifluxvariables.hh
@@ -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_;
};
diff --git a/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh b/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
index f7a22d5cfcdae5d43577022b5ad058b9d3c6ed42..c4a3450fe08ee6f00fbe5b68651d36c29a6ade58 100644
--- a/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnilocalresidual.hh
@@ -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_;
diff --git a/dumux/boxmodels/common/boxlocalresidual.hh b/dumux/boxmodels/common/boxlocalresidual.hh
index ff718f3aea10c6162ef65f98fe0e8e72fdb89c69..eaf68a38a3d1174ba2822ff580b8fd3275e013e1 100644
--- a/dumux/boxmodels/common/boxlocalresidual.hh
+++ b/dumux/boxmodels/common/boxlocalresidual.hh
@@ -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;
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
index 0563fd99828c84d793affa73a01ff137424ad47f..68ad0a23b5531d37fe56fd7fed9342db4f00c281 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
@@ -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);
}
diff --git a/dumux/material/spatialparameters/boxspatialparameters.hh b/dumux/material/spatialparameters/boxspatialparameters.hh
index 0573e80cc0156a921bea7ceb081d80846ddb0b0a..e1e9d44cbff854f279258e9781b2b513b735b537 100644
--- a/dumux/material/spatialparameters/boxspatialparameters.hh
+++ b/dumux/material/spatialparameters/boxspatialparameters.hh
@@ -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); }