diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 83b8b4c3261bb43035a6d3bea20c74c9083b5700..eb8b9577e0989f5cbf845d43e824aaeff605bbe7 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -317,6 +317,7 @@ $(schemes)/pointLinear/pointLinear.C $(schemes)/midPoint/midPoint.C $(schemes)/downwind/downwind.C $(schemes)/weighted/weighted.C +$(schemes)/weightedFlux/weightedFlux.C $(schemes)/cubic/cubic.C $(schemes)/skewCorrected/skewCorrectionVectors.C $(schemes)/skewCorrected/skewCorrected.C diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/weightedFlux/weightedFlux.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/weightedFlux/weightedFlux.C new file mode 100644 index 0000000000000000000000000000000000000000..05d53e1b90e9fc985ad5dfd4b34ac87d756c4ec7 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/weightedFlux/weightedFlux.C @@ -0,0 +1,239 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Norbert Weber, HZDR +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "weightedFlux.H" + +// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // + +template<class Type> +void Foam::weightedFlux<Type>::clearOut() +{ + deleteDemandDrivenData(oDelta_); + deleteDemandDrivenData(nDelta_); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template<class Type> +Foam::weightedFlux<Type>::~weightedFlux() +{ + clearOut(); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<class Type> +void Foam::weightedFlux<Type>::makeDeltas() const +{ + const fvMesh& mesh = this->mesh(); + + oDelta_ = new surfaceScalarField + ( + IOobject + ( + "oDelta", + mesh.pointsInstance(), + mesh + ), + mesh, + dimLength + ); + auto& oDelta = *oDelta_; + + nDelta_ = new surfaceScalarField + ( + IOobject + ( + "nDelta", + mesh.pointsInstance(), + mesh + ), + mesh, + dimLength + ); + auto& nDelta = *nDelta_; + + const labelUList& owner = mesh.owner(); + const labelUList& neighbour = mesh.neighbour(); + + const surfaceVectorField n = mesh.Sf()/mesh.magSf(); + + const vectorField& C = mesh.cellCentres(); + const vectorField& Cf = mesh.faceCentres(); + + // all distances are NORMAL to the face, + // as in the weighting factors in surfaceInterpolation.C + forAll(owner, facei) + { + oDelta[facei] = + mag(n[facei] & (C[owner[facei]] - Cf[facei])); + nDelta[facei] = + mag(n[facei] & (C[neighbour[facei]] - Cf[facei])); + } + + const fvPatchList& patches = mesh.boundary(); + + forAll(patches, patchi) + { + const fvPatch& currPatch = mesh.boundary()[patchi]; + + // Patch normal vector + const vectorField nPatch = currPatch.Sf()/currPatch.magSf(); + + // Processor patch + if (currPatch.coupled()) + { + const labelUList& pOwner = mesh.boundary()[patchi].faceCells(); + const vectorField& pCf = mesh.Cf().boundaryField()[patchi]; + + forAll(pOwner, facei) + { + const label own = pOwner[facei]; + + // All distances are NORMAL to the face + oDelta.boundaryFieldRef()[patchi][facei] = + mag(nPatch[facei] & (pCf[facei] - C[own])); + } + + // Weight = delta_neighbour / delta in ORTHOGONAL direction, + nDelta.boundaryFieldRef()[patchi] = + currPatch.weights()*oDelta.boundaryFieldRef()[patchi] + /(scalar(1) - currPatch.weights()); + } + else + { + const labelUList& pOwner = mesh.boundary()[patchi].faceCells(); + const vectorField& pCf = mesh.Cf().boundaryField()[patchi]; + + forAll(pOwner, facei) + { + const label own = pOwner[facei]; + + // All distances are NORMAL to the face! + oDelta.boundaryFieldRef()[patchi][facei] = + mag(nPatch[facei] & (pCf[facei] - C[own])); + + nDelta.boundaryFieldRef()[patchi][facei] = + mag(nPatch[facei] & (pCf[facei] - C[own])); + } + } + } +} + + +template<class Type> +Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>> +Foam::weightedFlux<Type>::interpolate +( + const GeometricField<Type, fvPatchField, volMesh>& vf +) const +{ + const fvMesh& mesh = vf.mesh(); + const surfaceScalarField& oDelta = weightedFlux<Type>::oDelta(); + const surfaceScalarField& nDelta = weightedFlux<Type>::nDelta(); + + auto tresult = tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>::New + ( + IOobject + ( + "weightedFlux::interpolate(" + vf.name() + ')', + mesh.time().timeName(), + mesh + ), + mesh, + vf.dimensions() + ); + auto& result = tresult.ref(); + + const labelUList& owner = mesh.owner(); + const labelUList& neighbour = mesh.neighbour(); + + forAll(result, facei) + { + const scalar sigmaDeltaO = sigma_[owner[facei]]/oDelta[facei]; + const scalar sigmaDeltaN = sigma_[neighbour[facei]]/nDelta[facei]; + + result[facei] = + (vf[owner[facei]]*sigmaDeltaO + vf[neighbour[facei]]*sigmaDeltaN) + /(sigmaDeltaO + sigmaDeltaN); + } + + // Interpolate patches + auto& bfld = result.boundaryFieldRef(); + + forAll(bfld, patchi) + { + fvsPatchField<Type>& pfld = bfld[patchi]; + + // If not coupled - simply copy the boundary values of the field + if (!pfld.coupled()) + { + pfld = vf.boundaryField()[patchi]; + } + else + { + // e.g. processor patches have to calculated separately + + const labelUList& pOwner = mesh.boundary()[patchi].faceCells(); + scalarField sigmaN = + sigma_.boundaryField()[patchi].patchNeighbourField(); + + Field<Type> vfO = vf.boundaryField()[patchi].patchInternalField(); + Field<Type> vfN = vf.boundaryField()[patchi].patchNeighbourField(); + + forAll(pOwner, facei) + { + const label own = pOwner[facei]; + + const scalar sigmaDeltaO = + sigma_[own]/oDelta.boundaryField()[patchi][facei]; + + const scalar sigmaDeltaN = + sigmaN[facei]/nDelta.boundaryField()[patchi][facei]; + + pfld[facei] = + (vfO[facei]*sigmaDeltaO + vfN[facei]*sigmaDeltaN) + /(sigmaDeltaO + sigmaDeltaN); + } + } + } + + return tresult; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makeSurfaceInterpolationScheme(weightedFlux) +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/weightedFlux/weightedFlux.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/weightedFlux/weightedFlux.H new file mode 100755 index 0000000000000000000000000000000000000000..4473c62e75aa67c224dc97c5a396dda90c3e6be2 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/weightedFlux/weightedFlux.H @@ -0,0 +1,218 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Norbert Weber, HZDR +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. + +Class + Foam::weightedFlux + +Description + Weighted flux interpolation scheme class. + + This scheme is used to compute fluxes with variable diffusivity or + conductivity, as e.g. + - a thermal flux: lambda*grad(T) + - a mass flux: D*grad(u) + - an electric current: -sigma*grad(potential) + + When using the Gauss theorem to compute a gradient, cell centred values + need to be interpolated to the faces. Using this scheme, temperature (T) + is weighted by thermal conductivity when being interpolated. Similarly, + velocity is weighted by diffusivity (D) and the electric potential by + the electric conductivity (sigma). Lambda, D or sigma are read from the + object registry - the names need to be specified in fvSchemes as e.g. + + \verbatim + gradSchemes + { + grad(T) Gauss weightedFlux "lambda"; + grad(u) Gauss weightedFlux "D"; + grad(potential) Gauss weightedFlux "sigma"; + } + \endverbatim + + For more details, see equation 16 and 17 in + \verbatim + Weber, N., Beckstein, P., Galindo, V., Starace, M. & Weier, T. (2018). + Electro-vortex flow simulation using coupled meshes. + Computers and Fluids 168, 101-109. + doi:10.1016/j.compfluid.2018.03.047 + https://arxiv.org/pdf/1707.06546.pdf + \endverbatim + +Note + For support, contact Norbert.Weber@hzdr.de + +SourceFiles + weightedFlux.C + +\*---------------------------------------------------------------------------*/ + +#ifndef weightedFlux_H +#define weightedFlux_H + +#include "surfaceInterpolationScheme.H" +#include "volFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class weightedFlux Declaration +\*---------------------------------------------------------------------------*/ + +template<class Type> +class weightedFlux +: + public surfaceInterpolationScheme<Type> +{ + // Private Data + + //- Const reference to step-wise pre-gradient factor field + const volScalarField& sigma_; + + // Demand-driven data + + //- Face to owner cell distance + mutable surfaceScalarField* oDelta_; + + //- Face to neighbour cell distance + mutable surfaceScalarField* nDelta_; + + + // Private Member Functions + + //- Compute face-owner and face-neighbour distance + void makeDeltas() const; + + //- No copy assignment + void operator=(const weightedFlux&) = delete; + + +protected: + + // Protected Member Functions + + // Storage management + + //- Clear all fields + void clearOut(); + + +public: + + //- Runtime type information + TypeName("weightedFlux"); + + + // Constructors + + //- Construct from Istream. + // The name of the flux field is read from the Istream and looked-up + // from the mesh objectRegistry + weightedFlux + ( + const fvMesh& mesh, + Istream& is + ) + : + surfaceInterpolationScheme<Type>(mesh), + sigma_(this->mesh().objectRegistry::template + lookupObject<volScalarField>(word(is))), + oDelta_(nullptr), + nDelta_(nullptr) + {} + + //- Construct from faceFlux and Istream + weightedFlux + ( + const fvMesh& mesh, + const surfaceScalarField& faceFlux, + Istream& is + ) + : + surfaceInterpolationScheme<Type>(mesh), + sigma_(this->mesh().objectRegistry::template + lookupObject<volScalarField>(word(is))), + oDelta_(nullptr), + nDelta_(nullptr) + {} + + + //- Destructor + ~weightedFlux(); + + + // Member Functions + + //- Return the interpolation weighting factors + tmp<surfaceScalarField> weights + ( + const GeometricField<Type, fvPatchField, volMesh>& + ) const + { + return this->mesh().surfaceInterpolation::weights(); + } + + //- Return the distance between face and owner cell + const surfaceScalarField& oDelta() const + { + if (!oDelta_) + { + makeDeltas(); + } + + return *oDelta_; + } + + //- Return the distance between face and neighbour cell + const surfaceScalarField& nDelta() const + { + if (!nDelta_) + { + makeDeltas(); + } + + return *nDelta_; + } + + //- Interpolate the cell values to faces + tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> + interpolate + ( + const GeometricField<Type, fvPatchField, volMesh>& vf + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //