From c1cf16b9946e26bc27931bbcd8e721901abb6530 Mon Sep 17 00:00:00 2001
From: Sina Ackermann <sina.ackermann@iws.uni-stuttgart.de>
Date: Mon, 6 Mar 2017 14:44:28 +0100
Subject: [PATCH] [staggeredGrid][freeflow] Adapt nonisothermal localresidual
and fluxvariables
* no test problem yet
---
dumux/freeflow/staggeredni/fluxvariables.hh | 114 +++++++++++---------
dumux/freeflow/staggeredni/localresidual.hh | 98 +----------------
2 files changed, 62 insertions(+), 150 deletions(-)
diff --git a/dumux/freeflow/staggeredni/fluxvariables.hh b/dumux/freeflow/staggeredni/fluxvariables.hh
index b8cfdf4ef4..d0acbf7f6f 100644
--- a/dumux/freeflow/staggeredni/fluxvariables.hh
+++ b/dumux/freeflow/staggeredni/fluxvariables.hh
@@ -65,13 +65,14 @@ class FreeFlowFluxVariablesImpl<TypeTag, false, true> : public FreeFlowFluxVaria
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); // TODO ?
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using IndexType = typename GridView::IndexSet::IndexType;
using Stencil = std::vector<IndexType>;
- using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
+// using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
// static constexpr auto numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
@@ -108,11 +109,11 @@ public:
const SubControlVolumeFace &scvf,
const FluxVariablesCache& fluxVarsCache)
{
-// CellCenterPrimaryVariables flux(0.0);
-//
-// flux += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, globalFaceVars, scvf);
-// flux += MolecularDiffusionType::diffusiveFluxForCellCenter(problem, fvGeometry, elemVolVars, scvf);
-// return flux;
+ CellCenterPrimaryVariables flux(0.0);
+
+ flux += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, globalFaceVars, scvf);
+ flux += diffusiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, scvf);
+ return flux;
}
private:
@@ -123,53 +124,60 @@ private:
const GlobalFaceVars& globalFaceVars,
const SubControlVolumeFace &scvf)
{
-// CellCenterPrimaryVariables flux(0.0);
-//
-// const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
-// const auto& insideVolVars = elemVolVars[insideScv];
-//
-// // if we are on an inflow/outflow boundary, use the volVars of the element itself
-// // TODO: catch neumann and outflow in localResidual's evalBoundary_()
-// bool isOutflow = false;
-// if(scvf.boundary())
-// {
-// const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
-// if(bcTypes.isOutflow(momentumBalanceIdx))
-// isOutflow = true;
-// }
-//
-// const auto& outsideVolVars = isOutflow ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
-//
-// const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
-//
-// const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity);
-// const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
-// const auto& downstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
-//
-// const Scalar upWindWeight = GET_PROP_VALUE(TypeTag, ImplicitUpwindWeight);
-// const Scalar upstreamDensity = useMoles ? upstreamVolVars.molarDensity() : upstreamVolVars.density();
-// const Scalar downstreamDensity = useMoles ? downstreamVolVars.molarDensity() : downstreamVolVars.density();
-//
-// for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-// {
-// // get equation index
-// const auto eqIdx = conti0EqIdx + compIdx;
-// if (eqIdx == replaceCompEqIdx)
-// continue;
-//
-// const Scalar upstreamFraction = useMoles ? upstreamVolVars.moleFraction(phaseIdx, compIdx) : upstreamVolVars.massFraction(phaseIdx, compIdx);
-// const Scalar downstreamFraction = useMoles ? downstreamVolVars.moleFraction(phaseIdx, compIdx) : downstreamVolVars.massFraction(phaseIdx, compIdx);
-//
-// flux[eqIdx] = (upWindWeight * upstreamDensity * upstreamFraction +
-// (1.0 - upWindWeight) * downstreamDensity * downstreamFraction)
-// * velocity;
-// }
-// // in case one balance is substituted by the total mass balance
-// if (replaceCompEqIdx < numComponents)
-// flux[replaceCompEqIdx] = (upWindWeight * upstreamDensity + (1.0 - upWindWeight) * downstreamDensity) * velocity;
-//
-// flux *= scvf.area() * sign(scvf.outerNormalScalar());
-// return flux;
+ CellCenterPrimaryVariables flux(0.0);
+
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& insideVolVars = elemVolVars[insideScv];
+
+ // if we are on an inflow/outflow boundary, use the volVars of the element itself
+ // TODO: catch neumann and outflow in localResidual's evalBoundary_()
+ bool isOutflow = false;
+ if(scvf.boundary())
+ {
+ const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+ if(bcTypes.isOutflow(momentumBalanceIdx))
+ isOutflow = true;
+ }
+
+ const auto& outsideVolVars = isOutflow ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
+
+ const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
+
+ const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity);
+ const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
+ const auto& downstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
+
+ const Scalar upWindWeight = GET_PROP_VALUE(TypeTag, ImplicitUpwindWeight);
+ const Scalar upstreamDensity = useMoles ? upstreamVolVars.molarDensity() : upstreamVolVars.density();
+ const Scalar downstreamDensity = useMoles ? downstreamVolVars.molarDensity() : downstreamVolVars.density();
+ const Scalar upstreamEnthalpy = upstreamVolVars.enthalpy();
+ const Scalar downstreamEnthalpy = downstreamVolVars.enthalpy();
+
+// flux[massBalanceIdx] = TODO??
+ flux[energyBalanceIdx] = (upWindWeight * upstreamDensity * upstreamEnthalpy
+ + (1.0 - upWindWeight) * downstreamDensity * downstreamEnthalpy)
+ * velocity;
+
+ flux *= scvf.area() * sign(scvf.outerNormalScalar());
+ return flux;
+ }
+
+ CellCenterPrimaryVariables diffusiveFluxForCellCenter_(const FluxVariables& fluxVars)
+ {
+ CellCenterPrimaryVariables flux(0.0);
+
+ // compute diffusive flux --> no diffusive flux (only 1 component)
+
+ // compute conductive flux
+ computeConductiveFlux_(flux, fluxVars);
+
+ return flux;
+ }
+
+ void computeConductiveFlux_(CellCenterPrimaryVariables& flux, FluxVariables& fluxVars)
+ {
+ flux[energyBalanceIdx] -= fluxVars.temperatureGrad() * fluxVars.face().normal
+ * (fluxVars.thermalConductivity() + fluxVars.thermalEddyConductivity());
}
};
diff --git a/dumux/freeflow/staggeredni/localresidual.hh b/dumux/freeflow/staggeredni/localresidual.hh
index 40d6ab4be4..d29b992386 100644
--- a/dumux/freeflow/staggeredni/localresidual.hh
+++ b/dumux/freeflow/staggeredni/localresidual.hh
@@ -78,10 +78,6 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false, true> : public Staggered
typename DofTypeIndices::FaceIdx faceIdx;
enum { // TODO adapt
- // grid and world dimension
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld,
-
pressureIdx = Indices::pressureIdx,
velocityIdx = Indices::velocityIdx,
@@ -110,105 +106,13 @@ public:
// const Scalar density = useMoles? volVars.molarDensity() : volVars.density();
// compute storage of mass
- storage[massBalanceIdx] = volVars.density(0);
+ storage[massBalanceIdx] = volVars.density(0); // TODO ParentType?
// compute the storage of energy
storage[energyBalanceIdx] = volVars.density(0) * volVars.internalEnergy(0);
return storage;
}
-
- // TODO implement advectiveFlux, conductiveFlux
-
- /*!
- * \brief Evaluates the convective energy flux
- * over a face of a sub-control volume and writes the result in
- * the flux vector. This method is called by computeFlux in the base class.
- *
- * \param flux The vector for the fluxes over the SCV/boundary face
- * \param fluxVars The flux variables at the current SCV/boundary face
- */
- void computeAdvectiveFlux(PrimaryVariables &flux,
- const FluxVariables &fluxVars) const
- {
- // call computation of the advective fluxes of the stokes model
- // (momentum and mass fluxes)
- ParentType::computeAdvectiveFlux(flux, fluxVars);
-
- // vertex data of the upstream and the downstream vertices
- const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx());
- const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx());
-
- Scalar tmp = fluxVars.normalVelocity();
- tmp *= (this->massUpwindWeight_ * up.density() * up.enthalpy()
- + (1.0 - this->massUpwindWeight_) * dn.density() * dn.enthalpy());
-
- flux[energyEqIdx] += tmp;
- Valgrind::CheckDefined(flux[energyEqIdx]);
- }
-
- /*!
- * \brief Evaluates the diffusive component energy flux
- * over the face of a sub-control volume.
- *
- * \param flux The vector for the fluxes over the SCV face
- * \param fluxVars The flux variables at the current SCV face
- */
- void computeDiffusiveFlux(PrimaryVariables &flux,
- const FluxVariables &fluxVars) const
- {
- // diffusive mass flux
- ParentType::computeDiffusiveFlux(flux, fluxVars);
-
- // conductive energy flux
- computeConductiveFlux(flux, fluxVars);
-
- // diffusive component energy flux
- Scalar sumDiffusiveFluxes = 0;
- for (int compIdx=0; compIdx<numComponents; compIdx++)
- {
- if (compIdx != phaseCompIdx)
- {
- Valgrind::CheckDefined(fluxVars.moleFractionGrad(compIdx));
- Valgrind::CheckDefined(fluxVars.face().normal);
- Valgrind::CheckDefined(fluxVars.diffusionCoeff(compIdx));
- Valgrind::CheckDefined(fluxVars.eddyDiffusivity());
- Valgrind::CheckDefined(fluxVars.molarDensity());
- Valgrind::CheckDefined(FluidSystem::molarMass(compIdx));
- Valgrind::CheckDefined(fluxVars.componentEnthalpy(compIdx));
- Scalar diffusiveFlux = fluxVars.moleFractionGrad(compIdx)
- * fluxVars.face().normal
- * (fluxVars.diffusionCoeff(compIdx) + fluxVars.eddyDiffusivity())
- * fluxVars.molarDensity();
- sumDiffusiveFluxes += diffusiveFlux;
- flux[energyEqIdx] -= diffusiveFlux * fluxVars.componentEnthalpy(compIdx)
- * FluidSystem::molarMass(compIdx); // Multiplied by molarMass [kg/mol] to convert from [mol/m^3 s] to [kg/m^3 s];
- }
- }
-
- // the diffusive flux of the phase component is the negative of the sum of the component fluxes
- flux[energyEqIdx] += sumDiffusiveFluxes * fluxVars.componentEnthalpy(phaseCompIdx)
- * FluidSystem::molarMass(phaseCompIdx); // Multiplied by molarMass [kg/mol] to convert from [mol/m^3 s] to [kg/m^3 s];
-
- Valgrind::CheckDefined(flux[energyEqIdx]);
- }
-
- /*!
- * \brief Evaluates the conductive energy flux
- * over the face of a sub-control volume.
- *
- * \param flux The vector for the fluxes over the SCV face
- * \param fluxVars The flux variables at the current SCV face
- */
- void computeConductiveFlux(PrimaryVariables &flux,
- const FluxVariables &fluxVars) const
- {
- // diffusive heat flux
- flux[energyEqIdx] -=
- fluxVars.temperatureGrad() * fluxVars.face().normal
- * (fluxVars.thermalConductivity() + fluxVars.thermalEddyConductivity());
- Valgrind::CheckDefined(flux[energyEqIdx]);
- }
};
} // end namespace
--
GitLab