exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Dynamic immersed boudnary with volume corrections, WIP

Browse source
This commit is contained in:
Hrvoje Jasak 2017-12-22 18:10:53 +00:00
parent f20d0995f9
commit 2ddecd498c
17 changed files with 810 additions and 269 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
src/immersedBoundary/immersedBoundary/Make/files
View file

@ -6,6 +6,8 @@ immersedBoundaryPointPatch/immersedBoundaryPointPatch.C
immersedBoundaryFvPatch/immersedBoundaryFvPatch.C
immersedBoundaryFvPatchField/immersedBoundaryFvPatchFields.C
immersedBoundaryFvsPatchField/immersedBoundaryFvsPatchFields.C
mixedIbFvPatchField/mixedIbFvPatchFields.C
movingImmersedBoundaryVelocity/movingImmersedBoundaryVelocityFvPatchVectorField.C
LIB = $(FOAM_LIBBIN)/libimmersedBoundary

132
src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C
View file

@ -35,7 +35,7 @@ License
namespace Foam
{
defineTypeNameAndDebug(immersedBoundaryFvPatch, 0);
defineTypeNameAndDebug(immersedBoundaryFvPatch, 1);
addToRunTimeSelectionTable(fvPatch, immersedBoundaryFvPatch, polyPatch);
}
@ -55,7 +55,6 @@ Foam::immersedBoundaryFvPatch::nonOrthogonalFactor_
void Foam::immersedBoundaryFvPatch::makeCf(slicedSurfaceVectorField& Cf) const
{
// Correct face centres by re-cutting and inserting the immersed patch
Info<< "Re-slicing Cf" << endl;
Cf.reset(ibPolyPatch_.correctedFaceCentres());
// Insert the patch data for the immersed boundary
@ -71,7 +70,6 @@ void Foam::immersedBoundaryFvPatch::makeCf(slicedSurfaceVectorField& Cf) const
void Foam::immersedBoundaryFvPatch::makeSf(slicedSurfaceVectorField& Sf) const
{
// Correct face centres by re-cutting and inserting the immersed patch
Info<< "Re-slicing Sf" << endl;
Sf.reset(ibPolyPatch_.correctedFaceAreas());
// Insert the patch data for the immersed boundary
@ -89,7 +87,6 @@ void Foam::immersedBoundaryFvPatch::makeSf(slicedSurfaceVectorField& Sf) const
void Foam::immersedBoundaryFvPatch::makeC(slicedVolVectorField& C) const
{
// Correct face centres by re-cutting and inserting the immersed patch
Info<< "Re-slicing C" << endl;
C.reset
(
ibPolyPatch_.correctedCellCentres(),
@ -109,11 +106,127 @@ void Foam::immersedBoundaryFvPatch::makeC(slicedVolVectorField& C) const
void Foam::immersedBoundaryFvPatch::makeV(scalarField& V) const
{
// Correct face centres by re-cutting and inserting the immersed patch
Info<< "Re-slicing V" << endl;
V = ibPolyPatch_.correctedCellVolumes();
}
void Foam::immersedBoundaryFvPatch::updatePhi
(
DimensionedField<scalar, volMesh>& V,
DimensionedField<scalar, volMesh>& V0,
surfaceScalarField& phi
) const
{
// Correct face fluxes for cut area and insert the immersed patch fluxes
const fvMesh& mesh = boundaryMesh().mesh();
scalar deltaT = mesh.time().deltaT().value();
scalar rDeltaT = 1.0/deltaT;
// Multiply the raw mesh motion flux with the masking function
scalarField sGamma =
mag(ibPolyPatch_.correctedFaceAreas())/mag(mesh.faceAreas());
// Scale internalField
phi.internalField() *= scalarField::subField(sGamma, mesh.nInternalFaces());
// Scale all other patches
forAll (mesh.boundary(), patchI)
{
if (!isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
{
phi.boundaryField()[patchI] *=
mesh.boundary()[patchI].patchSlice(sGamma);
}
}
// Immersed boundary patch
// Calculate the mesh motion flux from the old and new coordinate of
// triangular face centres and the time step dotted with the new face area
phi.boundaryField()[index()] =
(
ibPolyPatch_.motionDistance()
& ibPolyPatch_.correctedIbPatchFaceAreas()
)*rDeltaT;
// Check and adjust the immersed boundary space conservation law
// The mesh motion fluxes come from the actual mesh motion or the motion
// of the oimmersed boundary
// The new cell volumes come from the current mesh configuration
// Neither can be touched.
// The space conservation law will be satisfied by adjusting the old cell
// volume. HJ, 15/Dec/2017
// First sum up all the fluxes
scalarField divPhi(mesh.nCells(), 0);
const unallocLabelList& owner = mesh.owner();
const unallocLabelList& neighbour = mesh.neighbour();
const scalarField& phiIn = phi.internalField();
forAll(owner, faceI)
{
divPhi[owner[faceI]] += phiIn[faceI];
divPhi[neighbour[faceI]] -= phiIn[faceI];
}
// Add the mesh motion fluxes from all patches including immersed boundary
forAll(mesh.boundary(), patchI)
{
const unallocLabelList& pFaceCells =
mesh.boundary()[patchI].faceCells();
const scalarField& pssf = phi.boundaryField()[patchI];
forAll (pFaceCells, faceI)
{
divPhi[pFaceCells[faceI]] += pssf[faceI];
}
}
// Use corrected cell volume
scalarField& newVols = V.field();
scalarField& oldVols = V0.field();
// Multiply by the time-step size and add new volume
scalarField magDivPhi = mag((newVols - oldVols)*rDeltaT - divPhi);
// Note:
// The immersed boundary is now in the new position. Therefore, some
// cells that were cut are no longer in the contact with the IB, meaning
// that ALL cells need to be checked and corrected
// HJ, 22/Dec/2017
forAll (magDivPhi, cellI)
{
if (magDivPhi[cellI] > SMALL)
{
scalar corrOldVol = V[cellI] - divPhi[cellI]*deltaT;
// Info<< "Flux maneouvre for cell " << cellI << ": "
// << " V: " << V[cellI]
// << " V0: " << oldVols[cellI]
// << " divPhi: " << divPhi[cellI]
// << " Vcorr: " << corrOldVol
// << " error: " << magDivPhi[cellI]
// << " corr: "
// << (V[cellI] - Foam::max(corrOldVol, SMALL))*rDeltaT - divPhi[cellI]
// << endl;
if (corrOldVol < SMALL)
{
// Update new volume because old volume cannot carry
// the correction
newVols[cellI] = oldVols[cellI] + divPhi[cellI]*deltaT;
}
else
{
oldVols[cellI] = corrOldVol;
}
}
}
}
void Foam::immersedBoundaryFvPatch::makeCorrVecs(fvsPatchVectorField& cv) const
{
// Set patch non-orthogonality correction to zero
@ -122,7 +235,7 @@ void Foam::immersedBoundaryFvPatch::makeCorrVecs(fvsPatchVectorField& cv) const
vectorField& cvIn = const_cast<vectorField&>(cv.internalField());
const fvMesh& mesh = boundaryMesh().mesh();
// Get face addressing
const unallocLabelList& owner = mesh.owner();
const unallocLabelList& neighbour = mesh.neighbour();
@ -147,13 +260,6 @@ void Foam::immersedBoundaryFvPatch::makeCorrVecs(fvsPatchVectorField& cv) const
}
void Foam::immersedBoundaryFvPatch::movePoints()
{
Info<< "immersedBoundaryFvPatch::movePoints()" << endl;
ibPolyPatch_.clearOut();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::immersedBoundaryFvPatch::immersedBoundaryFvPatch

12
src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.H
View file

@ -106,14 +106,22 @@ protected:
//- Make mesh cell volumes
virtual void makeV(scalarField&) const;
//- Update mesh motion fluxes
virtual void updatePhi
(
DimensionedField<scalar, volMesh>& V,
DimensionedField<scalar, volMesh>& V0,
surfaceScalarField& phi
) const;
// Discretisation correction functions
//- Make patch face non-orthogonality correction vectors
virtual void makeCorrVecs(fvsPatchVectorField&) const;
//- Correct patches after moving points
virtual void movePoints();
// //- Correct patches after moving points
// virtual void movePoints();
public:

1
src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchField/immersedBoundaryFvPatchField.C
View file

@ -192,6 +192,7 @@ void immersedBoundaryFvPatchField<Type>::write(Ostream& os) const
// The value entry needs to be written with zero size
Field<Type>::null().writeEntry("value", os);
// this->writeEntry("value", os);
// Write VTK on master only
if (Pstream::master())

161
src/immersedBoundary/immersedBoundary/immersedBoundaryPolyPatch/immersedBoundaryPolyPatch.C
View file

@ -56,7 +56,7 @@ const Foam::debug::tolerancesSwitch
Foam::immersedBoundaryPolyPatch::liveFactor_
(
"immersedBoundaryLiveFactor",
1e-4
1e-3
);
@ -207,7 +207,7 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
// Used for sizing of dynamic list only
// HJ, 11/Dec/2017
label nIntersectedCells = 0;
// Go through the faces at the interface between a live and dead cell
// and mark the band of possible intersections
forAll (intersectedCell, cellI)
@ -412,7 +412,7 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
renumberedFace[fpI] = nIbPoints;
nIbPoints++;
}
// Record the face
unmergedFaces[nIbCells] = renumberedFace;
@ -460,7 +460,7 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
renumberedFace[fpI] = nIbPoints;
nIbPoints++;
}
// Record the face
unmergedFaces[nIbCells] = renumberedFace;
@ -505,7 +505,7 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
<< "Cannot find nearestTri for all points"
<< abort(FatalError);
}
// Build stand-alone patch
// Memory management
{
@ -541,18 +541,29 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
// Create IB patch from renumbered points and faces
ibPatchPtr_ = new standAlonePatch(ibPatchFaces, ibPatchPoints);
Info << "Writing immersed patch as VTK" << endl;
if (mesh.time().outputTime())
{
Info << "Writing immersed patch as VTK" << endl;
fileName fvPath(mesh.time().path()/"VTK");
mkDir(fvPath);
fileName fvPath(mesh.time().path()/"VTK");
mkDir(fvPath);
OStringStream outputFilename;
outputFilename << "immersed" << name();
fileName surfaceFileName
(
"immersed" + name() + "_live_"
+ Foam::name(boundaryMesh().mesh().time().timeIndex())
);
ibPatchPtr_->writeVTK(fvPath/fileName(outputFilename.str()));
ibPatchPtr_->writeVTK(fvPath/surfaceFileName);
outputFilename << "normals";
ibPatchPtr_->writeVTKNormals(fvPath/fileName(outputFilename.str()));
fileName normalsFileName
(
"normals" + name() + "_live_"
+ Foam::name(boundaryMesh().mesh().time().timeIndex())
);
ibPatchPtr_->writeVTKNormals(fvPath/normalsFileName);
}
}
// Count and collect dead cells
@ -969,12 +980,13 @@ void Foam::immersedBoundaryPolyPatch::calcCorrectedGeometry() const
&& mag(sumSf[ccc] + ibSf[cutCellI])/cutCellVolumes[cutCellI] > 1e-6
)
{
Info<< "Marooney Maneouvre for cell " << ccc
<< " error: " << mag(sumSf[ccc] + ibSf[cutCellI]) << " "
<< mag(sumSf[ccc] + ibSf[cutCellI])/cutCellVolumes[cutCellI]
<< " " << sumSf[ccc] << " "
<< " V: " << cutCellVolumes[cutCellI]
<< " Sf: " << ibSf[cutCellI] << endl;
// Info<< "Marooney Maneouvre for cell " << ccc
// << " error: " << mag(sumSf[ccc] + ibSf[cutCellI]) << " "
// << mag(sumSf[ccc] + ibSf[cutCellI])/cutCellVolumes[cutCellI]
// << " " << sumSf[ccc] << " "
// << " V: " << cutCellVolumes[cutCellI]
// << " Sf: " << ibSf[cutCellI] << endl;
ibSf[cutCellI] = -sumSf[ccc];
}
}
@ -991,14 +1003,20 @@ void Foam::immersedBoundaryPolyPatch::movePoints(const pointField& p)
(
"void immersedBoundaryPolyPatch::"
"movePoints(const pointField&) const"
) << "creating triSurface search algorithm"
) << "Moving mesh: immersedBoundary update"
<< endl;
}
// Handle motion of immersed boundary
clearOut();
// Handle motion of the mesh for new immersed boundary position
if (ibUpdateTimeIndex_ < boundaryMesh().mesh().time().timeIndex())
{
// New motion in the current time step. Clear
ibUpdateTimeIndex_ = boundaryMesh().mesh().time().timeIndex();
primitivePatch::movePoints(p);
clearOut();
}
polyPatch::movePoints(p);
}
@ -1043,7 +1061,8 @@ Foam::immersedBoundaryPolyPatch::immersedBoundaryPolyPatch
correctedFaceCentresPtr_(NULL),
correctedCellVolumesPtr_(NULL),
correctedFaceAreasPtr_(NULL),
correctedIbPatchFaceAreasPtr_(NULL)
correctedIbPatchFaceAreasPtr_(NULL),
oldIbPointsPtr_(NULL)
{}
@ -1087,7 +1106,8 @@ Foam::immersedBoundaryPolyPatch::immersedBoundaryPolyPatch
correctedFaceCentresPtr_(NULL),
correctedCellVolumesPtr_(NULL),
correctedFaceAreasPtr_(NULL),
correctedIbPatchFaceAreasPtr_(NULL)
correctedIbPatchFaceAreasPtr_(NULL),
oldIbPointsPtr_(NULL)
{
if (size() > 0)
{
@ -1145,7 +1165,8 @@ Foam::immersedBoundaryPolyPatch::immersedBoundaryPolyPatch
correctedFaceCentresPtr_(NULL),
correctedCellVolumesPtr_(NULL),
correctedFaceAreasPtr_(NULL),
correctedIbPatchFaceAreasPtr_(NULL)
correctedIbPatchFaceAreasPtr_(NULL),
oldIbPointsPtr_(NULL)
{}
@ -1188,7 +1209,8 @@ Foam::immersedBoundaryPolyPatch::immersedBoundaryPolyPatch
correctedFaceCentresPtr_(NULL),
correctedCellVolumesPtr_(NULL),
correctedFaceAreasPtr_(NULL),
correctedIbPatchFaceAreasPtr_(NULL)
correctedIbPatchFaceAreasPtr_(NULL),
oldIbPointsPtr_(NULL)
{}
@ -1197,6 +1219,8 @@ Foam::immersedBoundaryPolyPatch::immersedBoundaryPolyPatch
Foam::immersedBoundaryPolyPatch::~immersedBoundaryPolyPatch()
{
clearOut();
deleteDemandDrivenData(oldIbPointsPtr_);
}
@ -1393,6 +1417,47 @@ Foam::immersedBoundaryPolyPatch::correctedIbPatchFaceAreas() const
}
const Foam::pointField&
Foam::immersedBoundaryPolyPatch::oldIbPoints() const
{
if (!oldIbPointsPtr_)
{
// The mesh has never moved: old points are equal to current points
ibUpdateTimeIndex_ = boundaryMesh().mesh().time().timeIndex();
oldIbPointsPtr_ = new pointField(ibMesh_.points());
}
return *oldIbPointsPtr_;
}
Foam::tmp<Foam::vectorField>
Foam::immersedBoundaryPolyPatch::triMotionDistance() const
{
// Calculate the distance between new and old coordinates on
// the ibPatch face centres
// Calculate the motion on the triangular mesh face centres
return ibMesh_.coordinates()
- PrimitivePatch<labelledTri, List, const pointField&>
(
ibMesh_,
oldIbPoints()
).faceCentres();
}
Foam::tmp<Foam::vectorField>
Foam::immersedBoundaryPolyPatch::motionDistance() const
{
// Interpolate the values from tri surface using nearest triangle
return tmp<vectorField>
(
new vectorField(triMotionDistance(), nearestTri())
);
}
void Foam::immersedBoundaryPolyPatch::moveTriSurfacePoints
(
const pointField& p
@ -1418,26 +1483,42 @@ void Foam::immersedBoundaryPolyPatch::moveTriSurfacePoints
<< abort(FatalError);
}
if (ibUpdateTimeIndex_ < boundaryMesh().mesh().time().timeIndex())
{
// New motion in the current time step. Store old points
ibUpdateTimeIndex_ = boundaryMesh().mesh().time().timeIndex();
deleteDemandDrivenData(oldIbPointsPtr_);
Info<< "Storing old points for time index " << ibUpdateTimeIndex_
<< endl;
oldIbPointsPtr_ = new pointField(oldPoints);
}
Info<< "Moving immersed boundary points for patch " << name()
<< endl;
ibMesh_.movePoints(p);
fileName path(boundaryMesh().mesh().time().path()/"VTK");
if (boundaryMesh().mesh().time().outputTime())
{
fileName path(boundaryMesh().mesh().time().path()/"VTK");
mkDir(path);
ibMesh_.triSurface::write
(
path/
word
mkDir(path);
ibMesh_.triSurface::write
(
name() + "_"
+ Foam::name(boundaryMesh().mesh().time().timeIndex())
+ ".stl"
)
);
path/
word
(
name() + "_tri_"
+ Foam::name(boundaryMesh().mesh().time().timeIndex())
+ ".stl"
)
);
}
clearOut();
// Note: the IB patch is now in the new position, but the mesh has not
// been updated yet. movePoints() needs to be executed to update the
// fv mesh data
}
@ -1480,6 +1561,8 @@ void Foam::immersedBoundaryPolyPatch::clearOut() const
deleteDemandDrivenData(correctedCellVolumesPtr_);
deleteDemandDrivenData(correctedFaceAreasPtr_);
deleteDemandDrivenData(correctedIbPatchFaceAreasPtr_);
// Cannot delete old motion points. HJ, 10/Dec/2017
}

20
src/immersedBoundary/immersedBoundary/immersedBoundaryPolyPatch/immersedBoundaryPolyPatch.H
View file

@ -150,6 +150,12 @@ class immersedBoundaryPolyPatch
mutable vectorField* correctedIbPatchFaceAreasPtr_;
// Immersed boundary motion
//- Old points
mutable pointField* oldIbPointsPtr_;
// Private Member Functions
//- Check cell intersection
@ -350,11 +356,23 @@ public:
const vectorField& correctedIbPatchFaceAreas() const;
// Edit
// Immersed boundary motion
//- Correct patches after moving points
void moveTriSurfacePoints(const pointField& p);
//- Return old IB points
const pointField& oldIbPoints() const;
//- Motion distance vector for triangular immersed boundary faces
tmp<vectorField> triMotionDistance() const;
//- Motion distance vector for active immersed boundary faces
tmp<vectorField> motionDistance() const;
// Clear
//- Clear geometry
virtual void clearGeom();

52
src/immersedBoundary/immersedBoundary/include/createIbMasks.H
View file

@ -2,13 +2,18 @@
// createIbMask.H
// ~~~~~~~~~~~~
// Multiple IB patches currently not supported.
// This required IB geometry update to re-combine the corrected fields in
// Geometry correction functions in immersedBoundaryFvPatch
// HJ, 30/Nov/2017
Info<< "Create immersed boundary cell mask" << endl;
volScalarField cellIbMask
volScalarField gamma
(
IOobject
(
"cellIbMask",
"gamma",
runTime.timeName(),
mesh,
IOobject::NO_READ,
@ -19,11 +24,11 @@
);
Info<< "Create immersed boundary face mask" << endl;
surfaceScalarField faceIbMask
surfaceScalarField sGamma
(
IOobject
(
"faceIbMask",
"sGamma",
runTime.timeName(),
mesh,
IOobject::NO_READ,
@ -33,41 +38,4 @@
dimensionedScalar("one", dimless, 1)
);
// Multiple IB patches currently not supported.
// This required IB geometry update to re-combine the corrected fields in
// Geometry correction functions in immersedBoundaryFvPatch
// HJ, 30/Nov/2017
{
label nIbPatches = 0;
forAll (mesh.boundary(), patchI)
{
if (isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
{
Info<< "Found immersed boundary patch " << patchI
<< " named " << mesh.boundary()[patchI].name()
<< endl;
nIbPatches++;
const immersedBoundaryFvPatch& ibPatch =
refCast<const immersedBoundaryFvPatch>
(
mesh.boundary()[patchI]
);
cellIbMask = Foam::min(cellIbMask, ibPatch.gamma());
faceIbMask *= Foam::min(faceIbMask, ibPatch.sGamma());
}
}
if (nIbPatches > 1)
{
FatalErrorIn(args.executable())
<< "Multiple Immersed boundary patches currently not supported"
<< abort(FatalError);
}
}
// Evaluate boundary conditions for IB masks
cellIbMask.boundaryField().evaluateCoupled();
# include "updateIbMasks.H"

51
src/immersedBoundary/immersedBoundary/include/immersedBoundaryContinuityErrs.H
View file

@ -1,51 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Global
immersedBoundaryContinuityErrs
Description
Calculates and prints the continuity errors with
immersed boundary correction.
\*---------------------------------------------------------------------------*/
{
volScalarField contErr = fvc::div(faceIbMask*phi);
// volScalarField contErr = cellIbMask*fvc::div(phi);
sumLocalContErr = runTime.deltaT().value()*
mag(contErr)().weightedAverage(mesh.V()).value();
globalContErr = runTime.deltaT().value()*
contErr.weightedAverage(mesh.V()).value();
cumulativeContErr += globalContErr;
Info<< "IB time step continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr
<< endl;
}
// ************************************************************************* //

4
src/immersedBoundary/immersedBoundary/include/immersedBoundaryCourantNo.H
View file

@ -1,7 +1,7 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / O peration | Version: 3.2
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
@ -36,7 +36,7 @@ scalar velMag = 0.0;
if (mesh.nInternalFaces())
{
surfaceScalarField magPhi = mag(faceIbMask*phi);
surfaceScalarField magPhi = mag(sGamma*phi);
surfaceScalarField SfUfbyDelta =
mesh.surfaceInterpolation::deltaCoeffs()*magPhi;

78
src/immersedBoundary/immersedBoundary/include/updateIbMasks.H Normal file
View file

@ -0,0 +1,78 @@
{
// Multiple IB patches currently not supported.
// This required IB geometry update to re-combine the corrected fields in
// Geometry correction functions in immersedBoundaryFvPatch
// HJ, 30/Nov/2017
// Correct boundary patches to resize
gamma.correctBoundaryConditions();
sGamma.correctBoundaryConditions();
scalarField& gammaIn = gamma.internalField();
gammaIn = mesh.V()/mesh.cellVolumes();
scalarField& sGammaIn = sGamma.internalField();
const surfaceScalarField& magSf = mesh.magSf();
const scalarField magFaceAreas = mag(mesh.faceAreas());
sGammaIn = magSf.internalField()/
scalarField::subField(magFaceAreas, mesh.nInternalFaces());
forAll (mesh.boundary(), patchI)
{
if (!isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
{
sGamma.boundaryField()[patchI] =
magSf.boundaryField()[patchI]/
mesh.boundary()[patchI].patchSlice(magFaceAreas);
gamma.boundaryField()[patchI] =
sGamma.boundaryField()[patchI];
}
}
// Adjust the mask in dead faces and dead cells
forAll (mesh.boundary(), patchI)
{
if (isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
{
const immersedBoundaryFvPatch& ibPatch =
refCast<const immersedBoundaryFvPatch>
(
mesh.boundary()[patchI]
);
const labelList& deadCells = ibPatch.ibPolyPatch().deadCells();
forAll (deadCells, dcI)
{
gammaIn[deadCells[dcI]] = 0;
}
const labelList& deadFaces = ibPatch.ibPolyPatch().deadFaces();
forAll (deadFaces, dfI)
{
if (mesh.isInternalFace(deadFaces[dfI]))
{
sGammaIn[deadFaces[dfI]] = 0;
}
else
{
const label pI =
mesh.boundaryMesh().whichPatch(deadFaces[dfI]);
if (!isA<emptyFvPatch>(mesh.boundary()[pI]))
{
const label fI =
mesh.boundaryMesh()[pI].whichFace(deadFaces[dfI]);
sGamma.boundaryField()[pI][fI] = 0;
}
}
}
}
}
}

59
src/immersedBoundary/immersedBoundary/mixedIbFvPatchField/mixedIbFvPatchField.C
View file

@ -36,7 +36,7 @@ namespace Foam
template<class Type>
void mixedIbFvPatchField<Type>::updateIbValues()
{
// Interpolate the values form tri surface using nearest triangle
// Interpolate the values from tri surface using nearest triangle
const labelList& nt = ibPatch_.ibPolyPatch().nearestTri();
this->refValue() = Field<Type>(triValue_, nt);
@ -64,37 +64,6 @@ void mixedIbFvPatchField<Type>::setDeadValues()
}
template<class Type>
void mixedIbFvPatchField<Type>::correctDiag
(
fvMatrix<Type>& eqn
) const
{
scalarField& Diag = eqn.diag();
const labelList& deadCells = ibPatch_.ibPolyPatch().deadCells();
// Estimate diagonal in live cells
scalar liveDiag = 1;
if (deadCells.size() < Diag.size())
{
liveDiag = gSumMag(Diag)/(Diag.size() - deadCells.size());
// Correct for sign
liveDiag *= sign(gMax(Diag));
}
forAll (deadCells, cellI)
{
if (mag(Diag[deadCells[cellI]]) < SMALL)
{
Diag[deadCells[cellI]] = liveDiag;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
@ -289,31 +258,6 @@ void mixedIbFvPatchField<Type>::evaluate
}
// template<class Type>
// void mixedIbFvPatchField<Type>::manipulateMatrix
// (
// fvMatrix<Type>& eqn
// )
// {
// // Build matrix diagonal for cells where it is missing
// this->correctDiag(eqn);
// // Set values in IB cells
// Field<Type> polyPsi(eqn.psi(), ibPatch_.ibCells());
// eqn.setValues(ibPatch_.ibCells(), polyPsi);
// // Correct equation for dead cells
// Field<Type> deadCellsPsi
// (
// ibPatch_.deadCells().size(),
// deadValue_
// );
// eqn.setValues(ibPatch_.deadCells(), deadCellsPsi);
// fvPatchField<Type>::manipulateMatrix(eqn);
// }
template<class Type>
void mixedIbFvPatchField<Type>::write(Ostream& os) const
{
@ -331,6 +275,7 @@ void mixedIbFvPatchField<Type>::write(Ostream& os) const
// The value entry needs to be written with zero size
Field<Type>::null().writeEntry("value", os);
// this->writeEntry("value", os);
// Write VTK on master only
if (Pstream::master())

12
src/immersedBoundary/immersedBoundary/mixedIbFvPatchField/mixedIbFvPatchField.H
View file

@ -92,15 +92,6 @@ class mixedIbFvPatchField
void setDeadValues();
// Manipulate matrix
//- Check and correct zero diagonal in dead cells
void correctDiag
(
fvMatrix<Type>& eqn
) const;
public:
//- Runtime type information
@ -250,9 +241,6 @@ public:
const Pstream::commsTypes commsType = Pstream::blocking
);
// //- Manipulate matrix
// virtual void manipulateMatrix(fvMatrix<Type>& matrix);
// I-O

243
src/immersedBoundary/immersedBoundary/movingImmersedBoundaryVelocity/movingImmersedBoundaryVelocityFvPatchVectorField.C Normal file
View file

@ -0,0 +1,243 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "movingImmersedBoundaryVelocityFvPatchVectorField.H"
#include "addToRunTimeSelectionTable.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "fvcMeshPhi.H"
#include "surfaceWriter.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::movingImmersedBoundaryVelocityFvPatchVectorField::setDeadValues()
{
// Fix the value in dead cells
if (setDeadValue_)
{
const labelList& dc = ibPatch_.ibPolyPatch().deadCells();
// Get non-const access to internal field
vectorField& psiI = const_cast<vectorField&>(this->internalField());
forAll (dc, dcI)
{
psiI[dc[dcI]] = deadValue_;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::movingImmersedBoundaryVelocityFvPatchVectorField::
movingImmersedBoundaryVelocityFvPatchVectorField
(
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF
)
:
fixedValueFvPatchVectorField(p, iF),
ibPatch_(refCast<const immersedBoundaryFvPatch>(p)),
setDeadValue_(false),
deadValue_(vector::zero)
{}
Foam::movingImmersedBoundaryVelocityFvPatchVectorField::
movingImmersedBoundaryVelocityFvPatchVectorField
(
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF,
const dictionary& dict
)
:
fixedValueFvPatchVectorField(p, iF), // Do not read data
ibPatch_(refCast<const immersedBoundaryFvPatch>(p)),
setDeadValue_(dict.lookup("setDeadValue")),
deadValue_(dict.lookup("deadValue"))
{
// Evaluate with complete boundary velocity
vectorField::operator=
(
ibPatch_.ibPolyPatch().motionDistance()/
patch().boundaryMesh().mesh().time().deltaT().value()
);
}
Foam::movingImmersedBoundaryVelocityFvPatchVectorField::
movingImmersedBoundaryVelocityFvPatchVectorField
(
const movingImmersedBoundaryVelocityFvPatchVectorField& ptf,
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchVectorField(p, iF), // Do not data
ibPatch_(refCast<const immersedBoundaryFvPatch>(p)),
setDeadValue_(ptf.setDeadValue_),
deadValue_(ptf.deadValue_)
{
// Evaluate with complete boundary velocity
vectorField::operator=
(
ibPatch_.ibPolyPatch().motionDistance()/
patch().boundaryMesh().mesh().time().deltaT().value()
);
}
Foam::movingImmersedBoundaryVelocityFvPatchVectorField::
movingImmersedBoundaryVelocityFvPatchVectorField
(
const movingImmersedBoundaryVelocityFvPatchVectorField& ptf
)
:
fixedValueFvPatchVectorField(ptf),
ibPatch_(ptf.ibPatch()),
setDeadValue_(ptf.setDeadValue_),
deadValue_(ptf.deadValue_)
{}
Foam::movingImmersedBoundaryVelocityFvPatchVectorField::
movingImmersedBoundaryVelocityFvPatchVectorField
(
const movingImmersedBoundaryVelocityFvPatchVectorField& ptf,
const DimensionedField<vector, volMesh>& iF
)
:
fixedValueFvPatchVectorField(ptf, iF),
ibPatch_(ptf.ibPatch()),
setDeadValue_(ptf.setDeadValue_),
deadValue_(ptf.deadValue_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::movingImmersedBoundaryVelocityFvPatchVectorField::updateCoeffs()
{
if (updated())
{
return;
}
const fvMesh& mesh = dimensionedInternalField().mesh();
if (mesh.changing())
{
const fvPatch& p = patch();
// Get wall-parallel mesh motion velocity from immersed boundary
vectorField Up = ibPatch_.ibPolyPatch().motionDistance()/
mesh.time().deltaT().value();
const volVectorField& U =
mesh.lookupObject<volVectorField>
(
dimensionedInternalField().name()
);
scalarField phip =
p.patchField<surfaceScalarField, scalar>(fvc::meshPhi(U));
vectorField n = p.nf();
const scalarField& magSf = p.magSf();
scalarField Un = phip/(magSf + VSMALL);
// Adjust for surface-normal mesh motion flux
vectorField::operator=(Up + n*(Un - (n & Up)));
}
fixedValueFvPatchVectorField::updateCoeffs();
}
void Foam::movingImmersedBoundaryVelocityFvPatchVectorField::evaluate
(
const Pstream::commsTypes
)
{
// Set dead value
this->setDeadValues();
// Evaluate mixed condition
fixedValueFvPatchVectorField::evaluate();
}
void Foam::movingImmersedBoundaryVelocityFvPatchVectorField::write
(
Ostream& os
) const
{
fvPatchVectorField::write(os);
vectorField::null().writeEntry("value", os);
// writeEntry("value", os);
// Write VTK on master only
if (Pstream::master())
{
// Add parallel reduction of all faces and data to proc 0
// and write the whola patch together
// Write immersed boundary data as a vtk file
autoPtr<surfaceWriter<vector> > writerPtr =
surfaceWriter<vector>::New("vtk");
// Get the intersected patch
const standAlonePatch& ts = ibPatch_.ibPolyPatch().ibPatch();
writerPtr->write
(
this->dimensionedInternalField().path(),
ibPatch_.name(),
ts.points(),
ts,
this->dimensionedInternalField().name(),
*this,
surfaceWriterBase::FACE_DATA
);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
makePatchTypeField
(
fvPatchVectorField,
movingImmersedBoundaryVelocityFvPatchVectorField
);
} // End namespace Foam
// ************************************************************************* //

183
src/immersedBoundary/immersedBoundary/movingImmersedBoundaryVelocity/movingImmersedBoundaryVelocityFvPatchVectorField.H Normal file
View file

@ -0,0 +1,183 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::movingImmersedBoundaryVelocityFvPatchVectorField
Description
Foam::movingImmersedBoundaryVelocityFvPatchVectorField
SourceFiles
movingImmersedBoundaryVelocityFvPatchVectorField.C
\*---------------------------------------------------------------------------*/
#ifndef movingImmersedBoundaryVelocityFvPatchVectorField_H
#define movingImmersedBoundaryVelocityFvPatchVectorField_H
#include "fvPatchFields.H"
#include "fixedValueFvPatchFields.H"
#include "immersedBoundaryFvPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class movingImmersedBoundaryVelocityFvPatch Declaration
\*---------------------------------------------------------------------------*/
class movingImmersedBoundaryVelocityFvPatchVectorField
:
public fixedValueFvPatchVectorField
{
// Private data
//- Local reference to immersed boundary patch
const immersedBoundaryFvPatch& ibPatch_;
// Dead cell controls
//- Set dead cell value
Switch setDeadValue_;
//- Dead cell value
vector deadValue_;
// Private Member Functions
//- Set value in dead cells
void setDeadValues();
public:
//- Runtime type information
TypeName("movingImmersedBoundaryVelocity");
// Constructors
//- Construct from patch and internal field
movingImmersedBoundaryVelocityFvPatchVectorField
(
const fvPatch&,
const DimensionedField<vector, volMesh>&
);
//- Construct from patch, internal field and dictionary
movingImmersedBoundaryVelocityFvPatchVectorField
(
const fvPatch&,
const DimensionedField<vector, volMesh>&,
const dictionary&
);
//- Construct by mapping given
// movingImmersedBoundaryVelocityFvPatchVectorField
// onto a new patch
movingImmersedBoundaryVelocityFvPatchVectorField
(
const movingImmersedBoundaryVelocityFvPatchVectorField&,
const fvPatch&,
const DimensionedField<vector, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
movingImmersedBoundaryVelocityFvPatchVectorField
(
const movingImmersedBoundaryVelocityFvPatchVectorField&
);
//- Construct and return a clone
virtual tmp<fvPatchVectorField> clone() const
{
return tmp<fvPatchVectorField>
(
new movingImmersedBoundaryVelocityFvPatchVectorField(*this)
);
}
//- Construct as copy setting internal field reference
movingImmersedBoundaryVelocityFvPatchVectorField
(
const movingImmersedBoundaryVelocityFvPatchVectorField&,
const DimensionedField<vector, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchVectorField> clone
(
const DimensionedField<vector, volMesh>& iF
) const
{
return tmp<fvPatchVectorField>
(
new movingImmersedBoundaryVelocityFvPatchVectorField(*this, iF)
);
}
// Member functions
// Access
//- Return reference to immersed boundary patch
const immersedBoundaryFvPatch& ibPatch() const
{
return ibPatch_;
}
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Evaluate the patch field
virtual void evaluate
(
const Pstream::commsTypes commsType = Pstream::blocking
);
// I-O
//- Write
virtual void write(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

21
src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundarySolidBodyMotionFvMesh/immersedBoundarySolidBodyMotionFvMesh.C
View file

@ -41,9 +41,6 @@ namespace Foam
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::immersedBoundarySolidBodyMotionFvMesh::
@ -67,7 +64,8 @@ immersedBoundarySolidBodyMotionFvMesh
)
).subDict(typeName + "Coeffs")
),
ibMotions_()
ibMotions_(),
curTimeIndex_(-1)
{
// Read motion function for all regions
PtrList<entry> motionDicts(dynamicMeshCoeffs_.lookup("motionFunctions"));
@ -101,16 +99,29 @@ Foam::immersedBoundarySolidBodyMotionFvMesh::
bool Foam::immersedBoundarySolidBodyMotionFvMesh::update()
{
// Move points based on new motion
if (curTimeIndex_ < this->time().timeIndex())
{
// Grab old volumes before moving the mesh
setV0();
curTimeIndex_ = this->time().timeIndex();
}
forAll (ibMotions_, ibI)
{
ibMotions_[ibI].movePoints();
}
// Force quasi-topological cange to rebuild addressng on motion of the
// Force quasi-topological change to rebuild addressng on motion of the
// immersed boundary
// HJ, 12/Dec/2017
fvMesh::syncUpdateMesh();
// Execute dummy mesh motion for the background mesh
const pointField oldPoints = allPoints();
fvMesh::movePoints(oldPoints);
return true;
}

3
src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundarySolidBodyMotionFvMesh/immersedBoundarySolidBodyMotionFvMesh.H
View file

@ -61,6 +61,9 @@ class immersedBoundarySolidBodyMotionFvMesh
//- Immersed boundary motion control function
PtrList<movingImmersedBoundary> ibMotions_;
//- Current mesh motion time index
label curTimeIndex_;
// Private Member Functions

45
src/immersedBoundary/immersedBoundaryDynamicMesh/movingImmersedBoundary/movingImmersedBoundary.C
View file

@ -88,56 +88,11 @@ void Foam::movingImmersedBoundary::movePoints() const
immersedBoundaryPolyPatch& ibPatch =
const_cast<immersedBoundaryPolyPatch&>(cibPatch);
const vectorField oldIbPoints = ibPatch.ibMesh().coordinates();
// Move points
ibPatch.moveTriSurfacePoints
(
transform(sbmfPtr_->transformation(), refIbSurface_.points())
);
// Get non-const reference to velocity field
if (mesh().foundObject<volVectorField>("U"))
{
volVectorField& U = const_cast<volVectorField&>
(
mesh().lookupObject<volVectorField>("U")
);
if (isA<mixedIbFvPatchVectorField>(U.boundaryField()[patchID]))
{
// Get non-const reference to patch field
mixedIbFvPatchVectorField& ibPatchField =
refCast<mixedIbFvPatchVectorField>
(
U.boundaryField()[patchID]
);
// Set refValue_ to moving boundary velocity
ibPatchField.triValue() =
(ibPatch.ibMesh().coordinates() - oldIbPoints)/
mesh().time().deltaT().value();
}
else
{
InfoIn
(
"void movingImmersedBoundary::movePoints() const"
) << "Velocity field for patch "
<< name() << "is not a mixedIb type. "
<< "Immersed boundary velocity is not updated"
<< endl;
}
}
else
{
InfoIn
(
"void movingImmersedBoundary::movePoints() const"
) << "Velocity field not found for patch "
<< name() << ". Immersed boundary velocity is not updated"
<< endl;
}
}