interpolationCellPointWallModified: Restored interpolation method
This interpolation method was previously removed by commit fbf00209.
The intention of this method is to provide a slip-like wall boundary
condition for the velocity when interpolated to the location of a
Lagrangian element. This is difficult because any velocity which points
through the wall can cause a drag model and a rebound wall interaction
to "fight"; i.e., the drag pushes the particle to the wall, the wall
bounces it back. This can result in the program hanging.
This method extrapolates a vector field to the wall points and then
modifies the result so that it does not point through the wall. It does
this by rotating the vectors towards the (reversed) point normal. The
result is also scaled so that is reduced to zero if the necessary
rotation exceeds 90 degrees.
This provides an alternate resolution to bug report
https://bugs.openfoam.org/view.php?id=2826
This commit is contained in:
Will Bainbridge
@ -253,6 +253,7 @@ $(interpolation)/interpolationCellPatchConstrained/makeInterpolationCellPatchCon
|
||||
$(interpolation)/interpolationCellPoint/cellPointWeight/cellPointWeight.C
|
||||
$(interpolation)/interpolationCellPoint/makeInterpolationCellPoint.C
|
||||
$(interpolation)/interpolationCellPointFace/makeInterpolationCellPointFace.C
|
||||
$(interpolation)/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C
|
||||
$(interpolation)/interpolationPointMVC/pointMVCWeight.C
|
||||
$(interpolation)/interpolationPointMVC/makeInterpolationPointMVC.C
|
||||
|
||||
|
||||
@ -50,4 +50,19 @@ Foam::interpolationCellPoint<Type>::interpolationCellPoint
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
Foam::interpolationCellPoint<Type>::interpolationCellPoint
|
||||
(
|
||||
const GeometricField<Type, fvPatchField, volMesh>& psi,
|
||||
tmp<GeometricField<Type, pointPatchField, pointMesh>> psip
|
||||
)
|
||||
:
|
||||
interpolation<Type>(psi),
|
||||
psip_(psip)
|
||||
{
|
||||
// Uses cellPointWeight to do interpolation which needs tet decomposition
|
||||
(void)psi.mesh().tetBasePtIs();
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
@ -72,6 +72,13 @@ public:
|
||||
const GeometricField<Type, fvPatchField, volMesh>& psi
|
||||
);
|
||||
|
||||
//- Construct from components
|
||||
interpolationCellPoint
|
||||
(
|
||||
const GeometricField<Type, fvPatchField, volMesh>& psi,
|
||||
tmp<GeometricField<Type, pointPatchField, pointMesh>> psip
|
||||
);
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
|
||||
@ -0,0 +1,339 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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 "interpolationCellPointWallModified.H"
|
||||
#include "syncTools.H"
|
||||
#include "volPointInterpolation.H"
|
||||
#include "wallPolyPatch.H"
|
||||
#include "zeroGradientFvPatchFields.H"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
template<class Type>
|
||||
template<class TYPE>
|
||||
Foam::tmp<Foam::GeometricField<TYPE, Foam::pointPatchField, Foam::pointMesh>>
|
||||
Foam::interpolationCellPointWallModified<Type>::calcPointField
|
||||
(
|
||||
const GeometricField<TYPE, fvPatchField, volMesh>& psi
|
||||
) const
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< typeName << " interpolation is only defined for vector fields"
|
||||
<< exit(FatalError);
|
||||
|
||||
return tmp<GeometricField<TYPE, pointPatchField, pointMesh>>(nullptr);
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
Foam::tmp<Foam::pointVectorField>
|
||||
Foam::interpolationCellPointWallModified<Type>::calcPointField
|
||||
(
|
||||
const volVectorField& psi
|
||||
) const
|
||||
{
|
||||
using namespace constant::mathematical;
|
||||
|
||||
const fvMesh& mesh = psi.mesh();
|
||||
|
||||
|
||||
// Create the point field
|
||||
tmp<pointVectorField> tPsip
|
||||
(
|
||||
new pointVectorField
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"volPointInterpolateWallModified(" + psi.name() + ')',
|
||||
mesh.time().timeName(),
|
||||
mesh
|
||||
),
|
||||
pointMesh::New(mesh),
|
||||
dimensioned<Type>("zero", psi.dimensions(), Zero)
|
||||
)
|
||||
);
|
||||
pointVectorField& psip = tPsip.ref();
|
||||
|
||||
|
||||
// Interpolate to the points with wall patches extrapolated
|
||||
{
|
||||
wordList patchTypes(psi.boundaryField().size());
|
||||
forAll(patchTypes, patchi)
|
||||
{
|
||||
if (isA<wallPolyPatch>(mesh.boundaryMesh()[patchi]))
|
||||
{
|
||||
patchTypes[patchi] = zeroGradientFvPatchVectorField::typeName;
|
||||
}
|
||||
else
|
||||
{
|
||||
patchTypes[patchi] = calculatedFvPatchVectorField::typeName;
|
||||
}
|
||||
}
|
||||
volVectorField psiExtrapolated
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
psi.name() + "Extrapolated",
|
||||
mesh.time().timeName(),
|
||||
mesh
|
||||
),
|
||||
psi,
|
||||
patchTypes
|
||||
);
|
||||
psiExtrapolated.correctBoundaryConditions();
|
||||
volPointInterpolation::New(mesh).interpolate(psiExtrapolated, psip);
|
||||
}
|
||||
|
||||
|
||||
// Generate point normals across the entire boundary
|
||||
pointField pointNormals(mesh.nPoints(), vector::zero);
|
||||
{
|
||||
scalarField pointCount(mesh.nPoints(), 0);
|
||||
|
||||
forAll(mesh.boundaryMesh(), patchi)
|
||||
{
|
||||
const polyPatch& patch = mesh.boundaryMesh()[patchi];
|
||||
|
||||
forAll(patch, patchFacei)
|
||||
{
|
||||
const face& f = patch[patchFacei];
|
||||
const vector& n = patch.faceNormals()[patchFacei];
|
||||
|
||||
forAll(f, i)
|
||||
{
|
||||
pointNormals[f[i]] += n;
|
||||
pointCount[f[i]] += 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
syncTools::syncPointList
|
||||
(
|
||||
mesh,
|
||||
pointNormals,
|
||||
plusEqOp<vector>(),
|
||||
vector::zero
|
||||
);
|
||||
|
||||
syncTools::syncPointList
|
||||
(
|
||||
mesh,
|
||||
pointCount,
|
||||
plusEqOp<scalar>(),
|
||||
scalar(0)
|
||||
);
|
||||
|
||||
pointNormals /= max(pointCount, small);
|
||||
}
|
||||
|
||||
|
||||
// Calculate the rotation necessary from the vector to the (negative) point
|
||||
// normal needed to make the interpolated field point into the mesh
|
||||
scalarField theta0(mesh.nPoints(), -vGreat), theta1(mesh.nPoints(), vGreat);
|
||||
scalar maxVHatDotN = - vGreat;
|
||||
forAll(mesh.boundaryMesh(), patchi)
|
||||
{
|
||||
const polyPatch& patch = mesh.boundaryMesh()[patchi];
|
||||
if (isA<wallPolyPatch>(patch))
|
||||
{
|
||||
forAll(patch, patchFacei)
|
||||
{
|
||||
const label facei = patch.start() + patchFacei;
|
||||
const face& f = patch[patchFacei];
|
||||
|
||||
for (label i = 1; i < f.size() - 1; ++ i)
|
||||
{
|
||||
const triFace tri
|
||||
(
|
||||
tetIndices(mesh.faceOwner()[facei], facei, i)
|
||||
.faceTriIs(mesh)
|
||||
);
|
||||
|
||||
const vector n = tri.normal(mesh.points());
|
||||
|
||||
forAll(tri, triPointI)
|
||||
{
|
||||
const label pointi = tri[triPointI];
|
||||
|
||||
const vector& v = psip[pointi];
|
||||
|
||||
const scalar vHatDotN = normalised(v) & n;
|
||||
maxVHatDotN = max(maxVHatDotN, vHatDotN);
|
||||
|
||||
const vector a =
|
||||
normalised(v ^ pointNormals[pointi]);
|
||||
|
||||
const scalar C = v & n, S = (v ^ a) & n;
|
||||
|
||||
const scalar theta = atan2(C, - S);
|
||||
|
||||
theta0[pointi] = max(theta0[pointi], theta);
|
||||
theta1[pointi] = min(theta1[pointi], theta + pi);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (debug)
|
||||
{
|
||||
reduce(maxVHatDotN, maxOp<scalar>());
|
||||
Info<< typeName << ": Maximum in-to-wall dot product before = "
|
||||
<< maxVHatDotN << endl;
|
||||
}
|
||||
|
||||
syncTools::syncPointList
|
||||
(
|
||||
mesh,
|
||||
theta0,
|
||||
maxEqOp<scalar>(),
|
||||
scalar(0)
|
||||
);
|
||||
|
||||
syncTools::syncPointList
|
||||
(
|
||||
mesh,
|
||||
theta1,
|
||||
minEqOp<scalar>(),
|
||||
scalar(0)
|
||||
);
|
||||
|
||||
|
||||
if (debug > 1)
|
||||
{
|
||||
pointVectorField(psip.name() + "Before", psip).write();
|
||||
}
|
||||
|
||||
|
||||
// Apply the rotations so that the interpolated field points into the mesh
|
||||
forAll(mesh.boundaryMesh(), patchi)
|
||||
{
|
||||
const polyPatch& patch = mesh.boundaryMesh()[patchi];
|
||||
if (isA<wallPolyPatch>(patch))
|
||||
{
|
||||
forAll(patch.meshPoints(), patchPointi)
|
||||
{
|
||||
const label pointi = patch.meshPoints()[patchPointi];
|
||||
|
||||
vector& v = psip[pointi];
|
||||
|
||||
if (theta0[pointi] <= 0 && theta1[pointi] >= 0)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
if (theta0[pointi] >= theta1[pointi])
|
||||
{
|
||||
v = Zero;
|
||||
continue;
|
||||
}
|
||||
|
||||
const scalar theta =
|
||||
theta0[pointi] > 0 ? theta0[pointi] : theta1[pointi];
|
||||
|
||||
const scalar c = cos(theta), s = sin(theta);
|
||||
|
||||
const scalar scale = max(c, 0); // or mag(theta)/(pi/2) ...
|
||||
|
||||
const vector a = normalised(v ^ pointNormals[pointi]);
|
||||
|
||||
v = scale*(tensor::I*c - (*a)*s + sqr(a)*(1 - c)) & v;
|
||||
|
||||
theta0[pointi] = theta1[pointi] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Report the field-normal dot products
|
||||
if (debug)
|
||||
{
|
||||
maxVHatDotN = - vGreat;
|
||||
|
||||
forAll(mesh.boundaryMesh(), patchi)
|
||||
{
|
||||
const polyPatch& patch = mesh.boundaryMesh()[patchi];
|
||||
if (isA<wallPolyPatch>(patch))
|
||||
{
|
||||
forAll(patch, patchFacei)
|
||||
{
|
||||
const label facei = patch.start() + patchFacei;
|
||||
const face& f = patch[patchFacei];
|
||||
|
||||
for (label i = 1; i < f.size() - 1; ++ i)
|
||||
{
|
||||
const triFace tri
|
||||
(
|
||||
tetIndices(mesh.faceOwner()[facei], facei, i)
|
||||
.faceTriIs(mesh)
|
||||
);
|
||||
|
||||
const vector n = tri.normal(mesh.points());
|
||||
|
||||
forAll(tri, triPointI)
|
||||
{
|
||||
const label pointi = tri[triPointI];
|
||||
|
||||
const vector& v = psip[pointi];
|
||||
|
||||
const scalar vHatDotN = normalised(v) & n;
|
||||
maxVHatDotN = max(maxVHatDotN, vHatDotN);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
reduce(maxVHatDotN, maxOp<scalar>());
|
||||
Info<< typeName << ": Maximum in-to-wall dot product after = "
|
||||
<< maxVHatDotN << endl;
|
||||
}
|
||||
|
||||
|
||||
// Write the interpolated field
|
||||
if (debug > 1)
|
||||
{
|
||||
psip.write();
|
||||
}
|
||||
|
||||
|
||||
return tPsip;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
|
||||
|
||||
template<class Type>
|
||||
Foam::interpolationCellPointWallModified<Type>::
|
||||
interpolationCellPointWallModified
|
||||
(
|
||||
const GeometricField<Type, fvPatchField, volMesh>& psi
|
||||
)
|
||||
:
|
||||
interpolationCellPoint<Type>(psi, calcPointField(psi))
|
||||
{}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,106 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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::interpolationCellPointWallModified
|
||||
|
||||
Description
|
||||
As interpolationCellPoint, but with the point field modified on wall faces.
|
||||
|
||||
This method is defined only for vectors. The point field is extrapolated
|
||||
from the cells to the wall faces, and then rotated towards the (reverse)
|
||||
point normal so that the vectors do not point out of the domain. The result
|
||||
is also scaled so so if the necessary rotation exceeds 90 degrees, it is
|
||||
clamped to zero.
|
||||
|
||||
This prevents unresolvable drag-rebound couplings when applied to the
|
||||
velocity interpolation in a Lagrangian simulation.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef interpolationCellPointWallModified_H
|
||||
#define interpolationCellPointWallModified_H
|
||||
|
||||
#include "interpolationCellPoint.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class interpolationCellPointWallModified Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
template<class Type>
|
||||
class interpolationCellPointWallModified
|
||||
:
|
||||
public interpolationCellPoint<Type>
|
||||
{
|
||||
private:
|
||||
|
||||
// Private member functions
|
||||
|
||||
//- Compute the point field for a type. This just throws an error. This
|
||||
// interpolation method is only defined for vector fields.
|
||||
template <class TYPE>
|
||||
tmp<GeometricField<TYPE, pointPatchField, pointMesh>> calcPointField
|
||||
(
|
||||
const GeometricField<TYPE, fvPatchField, volMesh>& psi
|
||||
) const;
|
||||
|
||||
//- Compute the point field for a vector type
|
||||
tmp<pointVectorField> calcPointField(const volVectorField& psi) const;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
TypeName("cellPointWallModified");
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from components
|
||||
interpolationCellPointWallModified
|
||||
(
|
||||
const GeometricField<Type, fvPatchField, volMesh>& psi
|
||||
);
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#ifdef NoRepository
|
||||
#include "interpolationCellPointWallModified.C"
|
||||
#endif
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,35 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
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 "interpolationCellPointWallModified.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
makeInterpolation(interpolationCellPointWallModified);
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
Reference in New Issue
Block a user