Commit b601d7bf authored by Mark Olesen's avatar Mark Olesen
Browse files

Merge branch 'integration-feature-numerics' into 'develop'

Integration feature numerics

See merge request !308
parents 75ba4a07 3b6ba059
......@@ -14,23 +14,36 @@ volScalarField T
);
Info<< "Reading transportProperties\n" << endl;
Info<< "Reading diffusivity DT\n" << endl;
IOdictionary transportProperties
volScalarField DT
(
IOobject
(
"transportProperties",
runTime.constant(),
"DT",
runTime.timeName(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimViscosity, Zero)
);
Info<< "Reading diffusivity DT\n" << endl;
dimensionedScalar DT("DT", dimViscosity, transportProperties);
if (!DT.headerOk())
{
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
DT = dimensionedScalar("DT", dimViscosity, transportProperties);
}
#include "createFvOptions.H"
......@@ -41,6 +41,18 @@
gradT.component(vector::Z)
);
volVectorField DTgradT
(
IOobject
(
"flux",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
DT*gradT
);
runTime.write();
}
......@@ -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
......
/*---------------------------------------------------------------------------*\
========= |
\\ / 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)
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object DT;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -1 0 0 0 0];
internalField uniform 0;
boundaryField
{
"(top|bottom|sideWalls)"
{
type fixedGradient;
gradient uniform 0;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1912 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 1 0 0 0];
internalField uniform 273;
boundaryField