Robustness improvements and support for multiple immersed boundaries. Inno Gatin

This commit is contained in:
Hrvoje Jasak 2018-12-04 11:25:38 +00:00
parent 1a912d2914
commit 1b8d7cbce7
22 changed files with 2306 additions and 99 deletions

View file

@ -1,6 +1,7 @@
immersedPoly/immersedPoly.C
immersedPoly/distanceFunctions/triSurfaceDistance/triSurfaceDistance.C
immersedBoundaryPolyPatch/immersedBoundaryCorrectedMeshFields/immersedBoundaryCorrectedMeshFields.C
immersedBoundaryPolyPatch/immersedBoundaryPolyPatch.C
immersedBoundaryPointPatch/immersedBoundaryPointPatch.C
immersedBoundaryFvPatch/immersedBoundaryFvPatch.C

View file

@ -0,0 +1,174 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / 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 "immersedBoundaryFvcReconstruct.H"
#include "fvMesh.H"
#include "zeroGradientFvPatchFields.H"
#include "fv.H"
#include "surfaceInterpolate.H"
#include "fvcVolumeIntegrate.H"
#include "fvcSurfaceIntegrate.H"
#include "fvcAverage.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fvc
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
tmp
<
GeometricField
<
typename outerProduct<vector,Type>::type, fvPatchField, volMesh
>
>
reconstructIb
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
)
{
typedef typename outerProduct<vector, Type>::type GradType;
const fvMesh& mesh = ssf.mesh();
tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
(
new GeometricField<GradType, fvPatchField, volMesh>
(
IOobject
(
"volIntegrate(" + ssf.name() + ')',
ssf.instance(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
ssf.dimensions()/dimArea,
zeroGradientFvPatchField<GradType>::typeName
)
);
GeometricField<GradType, fvPatchField, volMesh>& reconField =
treconField();
// Note:
// 1) Reconstruction is only available in cell centres: there is no need
// to invert the tensor on the boundary
// 2) For boundaries, the only reconstructed data is the flux times
// normal. Based on this guess, boundary conditions can adjust
// patch values
// HJ, 12/Aug/2011
GeometricField<GradType, fvPatchField, volMesh> fluxTimesNormal =
surfaceSum((mesh.Sf()/mag(mesh.Sf()))*ssf);
// surfaceSum((mesh.Sf()/mesh.magSf())*ssf);
// Note: hinv inverse must be used to stabilise the inverse on bad meshes
// but it gives strange failures
// HJ, 19/Aug/2015
// Create the surface sum field
GeometricField<SymmTensor<scalar>, fvPatchField, volMesh> surfaceSumSqrSf =
surfaceSum(sqr(mesh.Sf())/mesh.magSf());
Field<SymmTensor<scalar> > surfaceSumSqrSfIn = surfaceSumSqrSf.internalField();
// Keep treck of the zero determinant cells
DynamicList<label> openCells(20);
// Check for zero determinant
forAll(surfaceSumSqrSfIn, cellI)
{
if (mag(det(surfaceSumSqrSfIn[cellI])) < VSMALL)
{
// Set the field to a unit tensor
surfaceSumSqrSfIn[cellI] = SymmTensor<scalar>::I;
// Store the cell ID for later correction of the field
openCells.append(cellI);
}
}
reconField.internalField() =
(
inv
(
surfaceSumSqrSfIn
)
& fluxTimesNormal.internalField()
);
// Set the zero determinant cells reconstructed field to zero
forAll(openCells, cellI)
{
reconField.internalField()[cellI] = GradType::zero;
}
reconField.boundaryField() = fluxTimesNormal.boundaryField();
treconField().correctBoundaryConditions();
return treconField;
}
template<class Type>
tmp
<
GeometricField
<
typename outerProduct<vector, Type>::type, fvPatchField, volMesh
>
>
reconstructIb
(
const tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >& tssf
)
{
typedef typename outerProduct<vector, Type>::type GradType;
tmp<GeometricField<GradType, fvPatchField, volMesh> > tvf
(
fvc::reconstructIb(tssf())
);
tssf.clear();
return tvf;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fvc
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,90 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / 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/>.
InNamespace
Foam::fvc
Description
Reconstruct volField from a face flux field.
SourceFiles
immersedBoundaryFvcReconstruct.C
\*---------------------------------------------------------------------------*/
#ifndef immersedBoundaryFvcReconstruct_H
#define immersedBoundaryFvcReconstruct_H
#include "volFieldsFwd.H"
#include "surfaceFieldsFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Namespace fvc functions Declaration
\*---------------------------------------------------------------------------*/
namespace fvc
{
template<class Type>
tmp
<
GeometricField
<typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
> reconstructIb
(
const GeometricField<Type, fvsPatchField, surfaceMesh>&
);
template<class Type>
tmp
<
GeometricField
<typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
> reconstructIb
(
const tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >&
);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "immersedBoundaryFvcReconstruct.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -60,10 +60,28 @@ void Foam::immersedBoundaryFvPatch::makeCf(slicedSurfaceVectorField& Cf) const
// Insert the patch data for the immersed boundary
// Note: use the face centres from the stand-alone patch within the IB
// HJ, 30/Nov/2017
Cf.boundaryField()[index()].UList::operator=
// NOTE: Loop though all ib boundaries in order to update the fields
// overwritten by reset. Needed for multiple IB patches (IG 2/Nov/2018)
const fvMesh& mesh = boundaryMesh().mesh();
forAll(mesh.boundary(), patchI)
{
if(isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
{
const immersedBoundaryFvPatch& curIbPatch = refCast
<
const immersedBoundaryFvPatch
>(mesh.boundary()[patchI]);
Cf.boundaryField()[patchI].UList::operator=
(
vectorField::subField(ibPolyPatch_.ibPatch().faceCentres(), size())
vectorField::subField
(
curIbPatch.ibPolyPatch().ibPatch().faceCentres(),
curIbPatch.size()
)
);
}
}
}
@ -76,10 +94,28 @@ void Foam::immersedBoundaryFvPatch::makeSf(slicedSurfaceVectorField& Sf) const
// Note: use the corrected face areas from immersed boundary instead of
// the stand-alone patch areas within the IB
// HJ, 30/Nov/2017
Sf.boundaryField()[index()].UList::operator=
// NOTE: Loop though all ib boundaries in order to update the fields
// overwritten by reset. Needed for multiple IB patches (IG 2/Nov/2018)
const fvMesh& mesh = boundaryMesh().mesh();
forAll(mesh.boundary(), patchI)
{
if(isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
{
const immersedBoundaryFvPatch& curIbPatch = refCast
<
const immersedBoundaryFvPatch
>(mesh.boundary()[patchI]);
Sf.boundaryField()[patchI].UList::operator=
(
vectorField::subField(ibPolyPatch_.correctedIbPatchFaceAreas(), size())
vectorField::subField
(
curIbPatch.ibPolyPatch().correctedIbPatchFaceAreas(),
curIbPatch.size()
)
);
}
}
}
@ -95,10 +131,28 @@ void Foam::immersedBoundaryFvPatch::makeC(slicedVolVectorField& C) const
// Insert the patch data for the immersed boundary
// Note: use the face centres from the stand-alone patch within the IB
// HJ, 30/Nov/2017
C.boundaryField()[index()].UList::operator=
// NOTE: Loop though all ib boundaries in order to update the fields
// overwritten by reset. Needed for multiple IB patches (IG 2/Nov/2018)
const fvMesh& mesh = boundaryMesh().mesh();
forAll(mesh.boundary(), patchI)
{
if(isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
{
const immersedBoundaryFvPatch& curIbPatch = refCast
<
const immersedBoundaryFvPatch
>(mesh.boundary()[patchI]);
C.boundaryField()[patchI].UList::operator=
(
vectorField::subField(ibPolyPatch_.ibPatch().faceCentres(), size())
vectorField::subField
(
curIbPatch.ibPolyPatch().ibPatch().faceCentres(),
curIbPatch.size()
)
);
}
}
}
@ -177,11 +231,16 @@ void Foam::immersedBoundaryFvPatch::updatePhi
const scalarField& pssf = phi.boundaryField()[patchI];
// Check for size since uninitialised ib patches can have zero size at
// this point (IG 7/Nov/2018)
if (pssf.size() > 0)
{
forAll (pFaceCells, faceI)
{
divPhi[pFaceCells[faceI]] += pssf[faceI];
}
}
}
// Use corrected cell volume
scalarField& newVols = V.field();
@ -246,7 +305,7 @@ void Foam::immersedBoundaryFvPatch::makeCorrVecs(fvsPatchVectorField& cv) const
if
(
gamma[owner[faceI]] < nonOrthogonalFactor_()
|| gamma[owner[faceI]] < nonOrthogonalFactor_()
|| gamma[neighbour[faceI]] < nonOrthogonalFactor_()
)
{
// Thin live cut. Reset correction vectors

View file

@ -40,6 +40,9 @@ Author
Zeljko Tukovic
Reorganisation and optimisation by Hrvoje Jasak
Contributors
Inno Gatin, FMENA, Zagreb.
SourceFiles
immersedBoundaryFvPatch.C
immersedBoundaryFvPatchTemplates.C

View file

@ -148,6 +148,7 @@ fixedValueIbFvPatchField<Type>::fixedValueIbFvPatchField
// Re-interpolate the data related to immersed boundary
this->updateIbValues();
this->setPatchType(ptf);
// On creation of the mapped field, the internal field is dummy and
// cannot be used. Initialise the value to avoid errors
@ -170,7 +171,9 @@ fixedValueIbFvPatchField<Type>::fixedValueIbFvPatchField
ptf.deadValue()
),
triValue_(ptf.triValue())
{}
{
this->setPatchType(ptf);
}
template<class Type>
@ -188,7 +191,9 @@ fixedValueIbFvPatchField<Type>::fixedValueIbFvPatchField
ptf.deadValue()
),
triValue_(ptf.triValue())
{}
{
this->setPatchType(ptf);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -258,6 +263,7 @@ void fixedValueIbFvPatchField<Type>::write(Ostream& os) const
// Resolve post-processing issues. HJ, 1/Dec/2017
fvPatchField<Type>::write(os);
triValue_.writeEntry("triValue", os);
immersedBoundaryFieldBase<Type>::writeDeadData(os);
// The value entry needs to be written with zero size
Field<Type>::null().writeEntry("value", os);

View file

@ -80,10 +80,10 @@ void Foam::immersedBoundaryFieldBase<Type>::writeField
) const
{
// Write VTK on master only
if (Pstream::master())
// if (Pstream::master())
{
// Add parallel reduction of all faces and data to proc 0
// and write the whole patch together
// and write the whole patch together. To Do: HJ, 4/Dec/2018
// Write immersed boundary data as a vtk file
autoPtr<surfaceWriter> writerPtr = surfaceWriter::New("vtk");

View file

@ -0,0 +1,134 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / 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 "immersedBoundaryCorrectedMeshFields.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(immersedBoundaryCorrectedMeshFields, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::immersedBoundaryCorrectedMeshFields::immersedBoundaryCorrectedMeshFields
(
const polyMesh& mesh
)
:
MeshObject<polyMesh, immersedBoundaryCorrectedMeshFields>(mesh),
correctedCellCentresPtr_(NULL),
correctedFaceCentresPtr_(NULL),
correctedCellVolumesPtr_(NULL),
correctedFaceAreasPtr_(NULL),
curTopoIndex_(-1)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::immersedBoundaryCorrectedMeshFields::~immersedBoundaryCorrectedMeshFields()
{
deleteDemandDrivenData(correctedCellCentresPtr_);
deleteDemandDrivenData(correctedFaceCentresPtr_);
deleteDemandDrivenData(correctedCellVolumesPtr_);
deleteDemandDrivenData(correctedFaceAreasPtr_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::immersedBoundaryCorrectedMeshFields::clearOut
(
const label topoChangeIndex
) const
{
// Check if the mesh changed since the last call and delete data if so.
// The fields must not get reset by every IB poly patch, this is why this
// check is necessary
if (curTopoIndex_ < topoChangeIndex)
{
deleteDemandDrivenData(correctedCellCentresPtr_);
deleteDemandDrivenData(correctedFaceCentresPtr_);
deleteDemandDrivenData(correctedCellVolumesPtr_);
deleteDemandDrivenData(correctedFaceAreasPtr_);
curTopoIndex_ = topoChangeIndex;
}
}
Foam::vectorField&
Foam::immersedBoundaryCorrectedMeshFields::correctedCellCentres() const
{
if (!correctedCellCentresPtr_)
{
correctedCellCentresPtr_ = new vectorField(mesh().cellCentres());
}
return *correctedCellCentresPtr_;
}
Foam::vectorField&
Foam::immersedBoundaryCorrectedMeshFields::correctedFaceCentres() const
{
if (!correctedFaceCentresPtr_)
{
correctedFaceCentresPtr_ = new vectorField(mesh().faceCentres());
}
return *correctedFaceCentresPtr_;
}
Foam::scalarField&
Foam::immersedBoundaryCorrectedMeshFields::correctedCellVolumes() const
{
if (!correctedCellVolumesPtr_)
{
correctedCellVolumesPtr_ = new scalarField(mesh().cellVolumes());
}
return *correctedCellVolumesPtr_;
}
Foam::vectorField&
Foam::immersedBoundaryCorrectedMeshFields::correctedFaceAreas() const
{
if (!correctedFaceAreasPtr_)
{
correctedFaceAreasPtr_ = new vectorField(mesh().faceAreas());
}
return *correctedFaceAreasPtr_;
}
// ************************************************************************* //

View file

@ -0,0 +1,133 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / 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::immersedBoundaryCorrectedMeshFields
Description
MeshObjet derived class for global corrected data in immersed boundary poly
patch.
SourceFiles
immersedBoundaryCorrectedMeshFields.C
Author
Inno Gatin, FMENA, Zagreb. All rights reserved.
\*---------------------------------------------------------------------------*/
#ifndef immersedBoundaryCorrectedMeshFields_H
#define immersedBoundaryCorrectedMeshFields_H
#include "MeshObject.H"
#include "polyMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class immersedBoundaryCorrectedMeshFields Declaration
\*---------------------------------------------------------------------------*/
class immersedBoundaryCorrectedMeshFields
:
public MeshObject<polyMesh, immersedBoundaryCorrectedMeshFields>
{
// Private data
//- Corrected cell centres
mutable vectorField* correctedCellCentresPtr_;
//- Corrected face centres
mutable vectorField* correctedFaceCentresPtr_;
//- Corrected cell volumes
mutable scalarField* correctedCellVolumesPtr_;
//- Corrected face areas
mutable vectorField* correctedFaceAreasPtr_;
//- Current topo change index
mutable label curTopoIndex_;
public:
// Declare name of the class and its debug switch
TypeName("immersedBoundaryCorrectedMeshFields");
// Constructors
//- Construct given an polyMesh
explicit immersedBoundaryCorrectedMeshFields(const polyMesh&);
// Destructor
virtual ~immersedBoundaryCorrectedMeshFields();
//- Clear data
void clearOut(const label topoChangeIndex) const;
// Access
// Corrected mesh geometry
//- Corrected cell centres
vectorField& correctedCellCentres() const;
//- Corrected face centres
vectorField& correctedFaceCentres() const;
//- Corrected cell volumes
scalarField& correctedCellVolumes() const;
//- Corrected face areas
vectorField& correctedFaceAreas() const;
// Virtual dummy functions: immersed boundary motion is handled elswere
virtual bool movePoints() const
{
return true;
}
virtual bool updateMesh(const mapPolyMesh&) const
{
return true;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -57,7 +57,7 @@ const Foam::debug::tolerancesSwitch
Foam::immersedBoundaryPolyPatch::liveFactor_
(
"immersedBoundaryLiveFactor",
1e-3
1e-6
);
@ -70,8 +70,9 @@ Foam::vector Foam::immersedBoundaryPolyPatch::cellSpan
{
const polyMesh& mesh = boundaryMesh().mesh();
// Calculate span from the bounding box size
const scalar delta = 3*cmptMax
// Calculate span from the bounding box size (prefactor is arbitrary, IG
// 10/Nov/2018)
const scalar delta = 20*cmptMax
(
boundBox
(
@ -245,12 +246,100 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
intersectedCell[cellI] = immersedPoly::DRY;
}
else if (foundInside && foundOutside)
{
// Get span
const vector span = cellSpan(cellI);
// If the nearest triangle cannot be found within span than this is
// most probably a tri surface search error. Mark unknown and check
// later. (IG 22/Nov/2018)
if (tss.nearest(C[cellI], span/20.0).index() == -1)
{
intersectedCell[cellI] = immersedPoly::UNKNOWN;
}
else
{
// Intersected cell
intersectedCell[cellI] = immersedPoly::CUT;
nIntersectedCells++;
}
}
}
// Do a check of the cells selected for cutting but not within the span of
// the tri surface. The cause of this can either be a stl that is not
// perfect or ther was an error in the inside/outside tri-search for other
// reasons. Look at the neigbours that are not CUT and assign their status.
const cellList& meshCells = mesh.cells();
forAll (intersectedCell, cellI)
{
if (intersectedCell[cellI] == immersedPoly::UNKNOWN)
{
// Check the neigbours
const cell& curCell = meshCells[cellI];
Switch foundWetNei = false;
Switch foundDryNei = false;
forAll (curCell, faceI)
{
// Only do the check for internal faces. If the face is boundary
// face then there is nothing to do. NOTE: parallelisation
// needed?
if (mesh.isInternalFace(curCell[faceI]))
{
label own = intersectedCell[owner[curCell[faceI]]];
label nei = intersectedCell[neighbour[curCell[faceI]]];
if
(
(nei == immersedPoly::DRY)
|| (own == immersedPoly::DRY)
)
{
foundDryNei = true;
}
if
(
(nei == immersedPoly::WET)
|| (own == immersedPoly::WET)
)
{
foundWetNei = true;
}
}
}
if (foundWetNei && !foundDryNei)
{
intersectedCell[cellI] = immersedPoly::WET;
}
else if (!foundWetNei && foundDryNei)
{
intersectedCell[cellI] = immersedPoly::DRY;
}
else
{
// There are either no wet or dry negbours or there are both.
// This should not be possible. NOTE: the check is not
// parallelised and this can theoretically lead to failiures in
// strange arrangaments.
// Issue a warning, mark CUT and hope for the best.
// (IG 22/Nov/2018)
WarningIn
(
"immersedBoundaryPolyPatch::calcImmersedBoundary() const"
)
<< "Cannot find wet or dry neigbours! Cell C:"
<< C[cellI]
<< " Neighbours: WET:" << foundWetNei
<< ", DRY:" << foundDryNei
<< endl;
intersectedCell[cellI] = immersedPoly::CUT;
nIntersectedCells++;
}
}
}
// Count all IB cells and faces for debug
labelList totalIbCount(4);
@ -312,7 +401,7 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
triSurfaceDistance dist
(
tss,
span,
2*span,
internalFlow(),
false // iterate intersection
);
@ -628,6 +717,7 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
// Record the nearest triangle to the face centre
nearestTri[nIbCells] =
tss.nearest(Cf[faceI], span).index();
nIbCells++;
}
}
@ -1085,11 +1175,7 @@ void Foam::immersedBoundaryPolyPatch::calcCorrectedGeometry() const
if
(
correctedCellCentresPtr_
|| correctedFaceCentresPtr_
|| correctedCellVolumesPtr_
|| correctedFaceAreasPtr_
|| correctedIbPatchFaceAreasPtr_
correctedIbPatchFaceAreasPtr_
)
{
FatalErrorIn
@ -1102,18 +1188,11 @@ void Foam::immersedBoundaryPolyPatch::calcCorrectedGeometry() const
// Get mesh reference
const polyMesh& mesh = boundaryMesh().mesh();
// Copy basic fields for corrections
correctedCellCentresPtr_ = new vectorField(mesh.cellCentres());
vectorField& C = *correctedCellCentresPtr_;
correctedFaceCentresPtr_ = new vectorField(mesh.faceCentres());
vectorField& Cf = *correctedFaceCentresPtr_;
correctedCellVolumesPtr_ = new scalarField(mesh.cellVolumes());
scalarField& V = *correctedCellVolumesPtr_;
correctedFaceAreasPtr_ = new vectorField(mesh.faceAreas());
vectorField& Sf = *correctedFaceAreasPtr_;
// Get mesh geometry references from the MeshObject
vectorField& C = correctedFields_.correctedCellCentres();;
vectorField& Cf = correctedFields_.correctedFaceCentres();
scalarField& V = correctedFields_.correctedCellVolumes();
vectorField& Sf = correctedFields_.correctedFaceAreas();
// Initialise IB patch face areas with the areas of the stand-alone patch
// They will be corrected using the Marooney Maneouvre
@ -1277,11 +1356,9 @@ Foam::immersedBoundaryPolyPatch::immersedBoundaryPolyPatch
nearestTriPtr_(NULL),
deadCellsPtr_(NULL),
deadFacesPtr_(NULL),
correctedCellCentresPtr_(NULL),
correctedFaceCentresPtr_(NULL),
correctedCellVolumesPtr_(NULL),
correctedFaceAreasPtr_(NULL),
correctedFields_(immersedBoundaryCorrectedMeshFields::New(bm.mesh())),
correctedIbPatchFaceAreasPtr_(NULL),
topoChangeIndex_(-1),
oldIbPointsPtr_(NULL)
{}
@ -1322,11 +1399,9 @@ Foam::immersedBoundaryPolyPatch::immersedBoundaryPolyPatch
nearestTriPtr_(NULL),
deadCellsPtr_(NULL),
deadFacesPtr_(NULL),
correctedCellCentresPtr_(NULL),
correctedFaceCentresPtr_(NULL),
correctedCellVolumesPtr_(NULL),
correctedFaceAreasPtr_(NULL),
correctedFields_(immersedBoundaryCorrectedMeshFields::New(bm.mesh())),
correctedIbPatchFaceAreasPtr_(NULL),
topoChangeIndex_(-1),
oldIbPointsPtr_(NULL)
{
if (size() > 0)
@ -1387,11 +1462,9 @@ Foam::immersedBoundaryPolyPatch::immersedBoundaryPolyPatch
nearestTriPtr_(NULL),
deadCellsPtr_(NULL),
deadFacesPtr_(NULL),
correctedCellCentresPtr_(NULL),
correctedFaceCentresPtr_(NULL),
correctedCellVolumesPtr_(NULL),
correctedFaceAreasPtr_(NULL),
correctedFields_(immersedBoundaryCorrectedMeshFields::New(bm.mesh())),
correctedIbPatchFaceAreasPtr_(NULL),
topoChangeIndex_(-1),
oldIbPointsPtr_(NULL)
{}
@ -1430,11 +1503,15 @@ Foam::immersedBoundaryPolyPatch::immersedBoundaryPolyPatch
nearestTriPtr_(NULL),
deadCellsPtr_(NULL),
deadFacesPtr_(NULL),
correctedCellCentresPtr_(NULL),
correctedFaceCentresPtr_(NULL),
correctedCellVolumesPtr_(NULL),
correctedFaceAreasPtr_(NULL),
correctedFields_
(
immersedBoundaryCorrectedMeshFields::New
(
pp.boundaryMesh().mesh()
)
),
correctedIbPatchFaceAreasPtr_(NULL),
topoChangeIndex_(-1),
oldIbPointsPtr_(NULL)
{}
@ -1474,11 +1551,9 @@ Foam::immersedBoundaryPolyPatch::immersedBoundaryPolyPatch
nearestTriPtr_(NULL),
deadCellsPtr_(NULL),
deadFacesPtr_(NULL),
correctedCellCentresPtr_(NULL),
correctedFaceCentresPtr_(NULL),
correctedCellVolumesPtr_(NULL),
correctedFaceAreasPtr_(NULL),
correctedFields_(immersedBoundaryCorrectedMeshFields::New(bm.mesh())),
correctedIbPatchFaceAreasPtr_(NULL),
topoChangeIndex_(-1),
oldIbPointsPtr_(NULL)
{}
@ -1629,48 +1704,48 @@ Foam::immersedBoundaryPolyPatch::deadFaces() const
const Foam::vectorField&
Foam::immersedBoundaryPolyPatch::correctedCellCentres() const
{
if (!correctedCellCentresPtr_)
if (!correctedIbPatchFaceAreasPtr_)
{
calcCorrectedGeometry();
}
return *correctedCellCentresPtr_;
return correctedFields_.correctedCellCentres();
}
const Foam::vectorField&
Foam::immersedBoundaryPolyPatch::correctedFaceCentres() const
{
if (!correctedFaceCentresPtr_)
if (!correctedIbPatchFaceAreasPtr_)
{
calcCorrectedGeometry();
}
return *correctedFaceCentresPtr_;
return correctedFields_.correctedFaceCentres();
}
const Foam::scalarField&
Foam::immersedBoundaryPolyPatch::correctedCellVolumes() const
{
if (!correctedCellVolumesPtr_)
if (!correctedIbPatchFaceAreasPtr_)
{
calcCorrectedGeometry();
}
return *correctedCellVolumesPtr_;
return correctedFields_.correctedCellVolumes();
}
const Foam::vectorField&
Foam::immersedBoundaryPolyPatch::correctedFaceAreas() const
{
if (!correctedFaceAreasPtr_)
if (!correctedIbPatchFaceAreasPtr_)
{
calcCorrectedGeometry();
}
return *correctedFaceAreasPtr_;
return correctedFields_.correctedFaceAreas();
}
@ -1825,12 +1900,14 @@ void Foam::immersedBoundaryPolyPatch::clearOut() const
deleteDemandDrivenData(deadCellsPtr_);
deleteDemandDrivenData(deadFacesPtr_);
deleteDemandDrivenData(correctedCellCentresPtr_);
deleteDemandDrivenData(correctedFaceCentresPtr_);
deleteDemandDrivenData(correctedCellVolumesPtr_);
deleteDemandDrivenData(correctedFaceAreasPtr_);
deleteDemandDrivenData(correctedIbPatchFaceAreasPtr_);
// Update topo change index
topoChangeIndex_++;
// Clear global data MeshObject
correctedFields_.clearOut(topoChangeIndex_);
// Cannot delete old motion points. HJ, 10/Dec/2017
}

View file

@ -32,6 +32,9 @@ Author
Reorganisation by Hrvoje Jasak
Complete rewrite of methodology, Hrvoje Jasak
Contributors
Inno Gatin, FMENA, Zagreb. All rights reserved.
SourceFiles
immersedBoundaryPolyPatch.C
@ -46,6 +49,8 @@ SourceFiles
#include "triSurfaceSearch.H"
#include "standAlonePatch.H"
#include "Switch.H"
#include "volFields.H"
#include "immersedBoundaryCorrectedMeshFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -134,21 +139,18 @@ class immersedBoundaryPolyPatch
// Corrected mesh geometry
//- Corrected cell centres
mutable vectorField* correctedCellCentresPtr_;
//- Corrected face centres
mutable vectorField* correctedFaceCentresPtr_;
//- Corrected cell volumes
mutable scalarField* correctedCellVolumesPtr_;
//- Corrected face areas
mutable vectorField* correctedFaceAreasPtr_;
// MeshObject derived object that holds corrected global field data
const immersedBoundaryCorrectedMeshFields& correctedFields_;
//- Corrected face areas for the immersed boundary patch
mutable vectorField* correctedIbPatchFaceAreasPtr_;
//- Topo change index, used to update the corrected fields
// MeshObject consistently.
// NOTE: time index cannot be used since multiple topo changes per
// time-step are possible (IG 28/OCT/2018)
mutable label topoChangeIndex_;
// Immersed boundary motion

View file

@ -353,6 +353,8 @@ Foam::face Foam::ImmersedCell<Distance>::createInternalFace
point dryPoint = vector::zero;
label nDry = 0;
label nUndecided = 0;
forAll (depth_, i)
{
if (depth_[i] > absTol_)
@ -365,6 +367,20 @@ Foam::face Foam::ImmersedCell<Distance>::createInternalFace
wetPoint += points_[i];
nWet++;
}
else
{
nUndecided++;
}
}
if(nUndecided == depth_.size())
{
FatalErrorIn
(
"ImmersedCell::createInternalFace(const label nIntersections)"
) << "All points lay on the tri surface, zero volume cell?"
<< nl << "Points: " << points_
<< abort(FatalError);
}
wetPoint /= nWet;
@ -372,7 +388,7 @@ Foam::face Foam::ImmersedCell<Distance>::createInternalFace
// Good direction points out of the wet cell
vector dir = dryPoint - wetPoint;
dir /= mag(dir);
dir /= mag(dir) + SMALL;
vector n = orderedInternalFace.normal(points_);
n /= mag(n);

View file

@ -20,6 +20,19 @@
sGammaIn = magSf.internalField()/
scalarField::subField(magFaceAreas, mesh.nInternalFaces());
const scalar maxGamma = max(gammaIn);
const scalar maxSGamma = min(sGammaIn);
if (maxGamma > 1 || maxSGamma > 1)
{
WarningIn
(
"updateIbMasks.H"
)
<< "Unbounded gamma: max(gamma) = "
<< maxGamma << ", max(sGamma) = " <<
<< maxSGamma << endl;
}
forAll (mesh.boundary(), patchI)
{
if (!isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))

View file

@ -3,4 +3,7 @@ wallFunctions/immersedBoundaryNutWallFunction/immersedBoundaryNutWallFunctionFvP
wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C
wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.C
turbulenceModels/incompressible/RAS/backwardsCompability/wallFunctions/backwardsCompatibilityIbWallFunctions.C
turbulenceModels/incompressible/RAS/immersedBoundaryKOmegaSST/immersedBoundaryKOmegaSST.C
LIB = $(FOAM_LIBBIN)/libimmersedBoundaryTurbulence

View file

@ -0,0 +1,336 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / 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 "backwardsCompatibilityIbWallFunctions.H"
#include "calculatedFvPatchField.H"
#include "nutkWallFunctionFvPatchScalarField.H"
#include "nutLowReWallFunctionFvPatchScalarField.H"
#include "epsilonWallFunctionFvPatchScalarField.H"
#include "kqRWallFunctionFvPatchField.H"
#include "RWallFunctionFvPatchSymmTensorField.H"
#include "omegaWallFunctionFvPatchScalarField.H"
#include "immersedBoundaryEpsilonWallFunctionFvPatchScalarField.H"
#include "immersedBoundaryKqRWallFunctionFvPatchField.H"
#include "immersedBoundaryOmegaWallFunctionFvPatchScalarField.H"
#include "immersedBoundaryNutWallFunctionFvPatchScalarField.H"
#include "immersedBoundaryFvPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
tmp<volScalarField> autoCreateIbNut
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
)
{
IOobject nutHeader
(
fieldName,
mesh.time().timeName(),
obj,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
);
if (nutHeader.headerOk())
{
return tmp<volScalarField>(new volScalarField(nutHeader, mesh));
}
else
{
Info<< "--> Creating " << fieldName
<< " to employ run-time selectable wall functions" << endl;
const fvBoundaryMesh& bm = mesh.boundary();
wordList nutBoundaryTypes(bm.size());
forAll(bm, patchI)
{
if (isA<immersedBoundaryFvPatch>(bm[patchI]) && bm[patchI].isWall())
{
nutBoundaryTypes[patchI] =
RASModels::immersedBoundaryNutWallFunctionFvPatchScalarField::typeName;
}
else if (bm[patchI].isWall())
{
nutBoundaryTypes[patchI] =
RASModels::nutkWallFunctionFvPatchScalarField::typeName;
}
else
{
nutBoundaryTypes[patchI] =
calculatedFvPatchField<scalar>::typeName;
}
}
tmp<volScalarField> nut
(
new volScalarField
(
IOobject
(
fieldName,
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar("zero", dimArea/dimTime, 0.0),
nutBoundaryTypes
)
);
Info<< " Writing new " << fieldName << endl;
nut().write();
return nut;
}
}
tmp<volScalarField> autoCreateIbNut
(
const word& fieldName,
const fvMesh& mesh
)
{
return autoCreateIbNut(fieldName, mesh, mesh);
}
tmp<volScalarField> autoCreateIbEpsilon
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
)
{
return
autoCreateWallFunctionField
<
scalar,
RASModels::epsilonWallFunctionFvPatchScalarField,
RASModels::immersedBoundaryEpsilonWallFunctionFvPatchScalarField
>
(
fieldName,
mesh,
obj
);
}
tmp<volScalarField> autoCreateIbEpsilon
(
const word& fieldName,
const fvMesh& mesh
)
{
return
autoCreateWallFunctionField
<
scalar,
RASModels::epsilonWallFunctionFvPatchScalarField,
RASModels::immersedBoundaryEpsilonWallFunctionFvPatchScalarField
>
(
fieldName,
mesh,
mesh
);
}
tmp<volScalarField> autoCreateIbOmega
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
)
{
return
autoCreateWallFunctionField
<
scalar,
RASModels::omegaWallFunctionFvPatchScalarField,
RASModels::immersedBoundaryOmegaWallFunctionFvPatchScalarField
>
(
fieldName,
mesh,
obj
);
}
tmp<volScalarField> autoCreateIbOmega
(
const word& fieldName,
const fvMesh& mesh
)
{
return
autoCreateWallFunctionField
<
scalar,
RASModels::omegaWallFunctionFvPatchScalarField,
RASModels::immersedBoundaryOmegaWallFunctionFvPatchScalarField
>
(
fieldName,
mesh,
mesh
);
}
tmp<volScalarField> autoCreateIbK
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
)
{
return
autoCreateWallFunctionField
<
scalar,
RASModels::kqRWallFunctionFvPatchField<scalar>,
RASModels::immersedBoundaryKqRWallFunctionFvPatchField<scalar>
>
(
fieldName,
mesh,
obj
);
}
tmp<volScalarField> autoCreateIbK
(
const word& fieldName,
const fvMesh& mesh
)
{
return
autoCreateWallFunctionField
<
scalar,
RASModels::kqRWallFunctionFvPatchField<scalar>,
RASModels::immersedBoundaryKqRWallFunctionFvPatchField<scalar>
>
(
fieldName,
mesh,
mesh
);
}
tmp<volScalarField> autoCreateIbQ
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
)
{
return
autoCreateWallFunctionField
<
scalar,
RASModels::kqRWallFunctionFvPatchField<scalar>,
RASModels::immersedBoundaryKqRWallFunctionFvPatchField<scalar>
>
(
fieldName,
mesh,
obj
);
}
tmp<volScalarField> autoCreateIbQ
(
const word& fieldName,
const fvMesh& mesh
)
{
return
autoCreateWallFunctionField
<
scalar,
RASModels::kqRWallFunctionFvPatchField<scalar>,
RASModels::immersedBoundaryKqRWallFunctionFvPatchField<scalar>
>
(
fieldName,
mesh,
mesh
);
}
tmp<volSymmTensorField> autoCreateIbR
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
)
{
return
autoCreateWallFunctionField
<
symmTensor,
RASModels::kqRWallFunctionFvPatchField<symmTensor>,
RASModels::immersedBoundaryKqRWallFunctionFvPatchField<symmTensor>
>
(
fieldName,
mesh,
obj
);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace incompressible
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,158 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / 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::incompressible
Description
Auto creation of fields to provide backwards compatibility with
runtime selectable wall functions
SourceFiles
backwardsCompatibilityIbWallFunctions.C
backwardsCompatibilityIbWallFunctionsTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef backwardsCompatibilityIbWallFunctions_H
#define backwardsCompatibilityIbWallFunctions_H
#include "fvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
//- nut
tmp<volScalarField> autoCreateIbNut
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
);
tmp<volScalarField> autoCreateIbNut
(
const word& fieldName,
const fvMesh& mesh
);
//- epsilon
tmp<volScalarField> autoCreateIbEpsilon
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
);
tmp<volScalarField> autoCreateIbEpsilon
(
const word& fieldName,
const fvMesh& mesh
);
//- omega
tmp<volScalarField> autoCreateIbOmega
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
);
tmp<volScalarField> autoCreateIbOmega
(
const word& fieldName,
const fvMesh& mesh
);
//- k
tmp<volScalarField> autoCreateIbK
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
);
tmp<volScalarField> autoCreateIbK
(
const word& fieldName,
const fvMesh& mesh
);
//- Q
tmp<volScalarField> autoCreateIbQ
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
);
tmp<volScalarField> autoCreateIbQ
(
const word& fieldName,
const fvMesh& mesh
);
//- R
tmp<volSymmTensorField> autoCreateIbR
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
);
tmp<volSymmTensorField> autoCreateIbR
(
const word& fieldName,
const fvMesh& mesh
);
//- Helper function to create the new field
template<class Type, class PatchType, class ibPatchName >
tmp<GeometricField<Type, fvPatchField, volMesh> >
autoCreateIbWallFunctionField
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace incompressible
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "backwardsCompatibilityIbWallFunctionsTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -0,0 +1,184 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / 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 "backwardsCompatibilityIbWallFunctions.H"
#include "foamTime.H"
#include "OSspecific.H"
#include "immersedBoundaryFvPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class PatchType, class ibPatchType>
tmp<GeometricField<Type, fvPatchField, volMesh> >
autoCreateWallFunctionField
(
const word& fieldName,
const fvMesh& mesh,
const objectRegistry& obj
)
{
IOobject nutHeader
(
"nut",
mesh.time().timeName(),
obj,
IOobject::MUST_READ
);
typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
if (nutHeader.headerOk())
{
return tmp<fieldType>
(
new fieldType
(
IOobject
(
fieldName,
mesh.time().timeName(),
obj,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
),
mesh
)
);
}
else
{
Info<< "--> Upgrading " << fieldName
<< " to employ run-time selectable wall functions" << endl;
// Read existing field
IOobject ioObj
(
fieldName,
mesh.time().timeName(),
obj,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
);
tmp<fieldType> fieldOrig
(
new fieldType
(
ioObj,
mesh
)
);
// rename file
Info<< " Backup original " << fieldName << " to "
<< fieldName << ".old" << endl;
mvBak(ioObj.objectPath(), "old");
PtrList<fvPatchField<Type> > newPatchFields(mesh.boundary().size());
const fvBoundaryMesh& bm = mesh.boundary();
forAll(newPatchFields, patchI)
{
if (isA<immersedBoundaryFvPatch>(bm[patchI]) && bm[patchI].isWall())
{
newPatchFields.set
(
patchI,
new ibPatchType
(
mesh.boundary()[patchI],
fieldOrig().dimensionedInternalField()
)
);
newPatchFields[patchI] == fieldOrig().boundaryField()[patchI];
}
else if (bm[patchI].isWall())
{
newPatchFields.set
(
patchI,
new PatchType
(
mesh.boundary()[patchI],
fieldOrig().dimensionedInternalField()
)
);
newPatchFields[patchI] == fieldOrig().boundaryField()[patchI];
}
else
{
newPatchFields.set
(
patchI,
fieldOrig().boundaryField()[patchI].clone()
);
}
}
tmp<fieldType> fieldNew
(
new fieldType
(
IOobject
(
fieldName,
mesh.time().timeName(),
obj,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
fieldOrig().dimensions(),
fieldOrig().internalField(),
newPatchFields
)
);
Info<< " Writing updated " << fieldName << endl;
fieldNew().write();
return fieldNew;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace incompressible
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,483 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / 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 "immersedBoundaryKOmegaSST.H"
#include "addToRunTimeSelectionTable.H"
#include "backwardsCompatibilityIbWallFunctions.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace RASModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(immersedBoundaryKOmegaSST, 0);
addToRunTimeSelectionTable(RASModel, immersedBoundaryKOmegaSST, dictionary);
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
tmp<volScalarField> immersedBoundaryKOmegaSST::F1
(
const volScalarField& CDkOmega
) const
{
volScalarField CDkOmegaPlus = max
(
CDkOmega,
dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
);
volScalarField arg1 = min
(
min
(
max
(
(scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
scalar(500)*nu()/(sqr(y_)*omega_)
),
(4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
),
scalar(10)
);
return tanh(pow4(arg1));
}
tmp<volScalarField> immersedBoundaryKOmegaSST::F2() const
{
volScalarField arg2 = min
(
max
(
(scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
scalar(500)*nu()/(sqr(y_)*omega_)
),
scalar(100)
);
return tanh(sqr(arg2));
}
tmp<volScalarField> immersedBoundaryKOmegaSST::F3() const
{
tmp<volScalarField> arg3 = min
(
150*nu()/(omega_*sqr(y_)),
scalar(10)
);
return 1 - tanh(pow4(arg3));
}
tmp<volScalarField> immersedBoundaryKOmegaSST::F23() const
{
tmp<volScalarField> f23(F2());
if (F3_)
{
f23() *= F3();
}
return f23;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
immersedBoundaryKOmegaSST::immersedBoundaryKOmegaSST
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName,
const word& modelName
)
:
RASModel(modelName, U, phi, transport, turbulenceModelName),
alphaK1_
(
dimensionedScalar::lookupOrAddToDict
(
"alphaK1",
coeffDict_,
0.85
)
),
alphaK2_
(
dimensionedScalar::lookupOrAddToDict
(
"alphaK2",
coeffDict_,
1.0
)
),
alphaOmega1_
(
dimensionedScalar::lookupOrAddToDict
(
"alphaOmega1",
coeffDict_,
0.5
)
),
alphaOmega2_
(
dimensionedScalar::lookupOrAddToDict
(
"alphaOmega2",
coeffDict_,
0.856
)
),
gamma1_
(
dimensionedScalar::lookupOrAddToDict
(
"gamma1",
coeffDict_,
5.0/9.0
)
),
gamma2_
(
dimensionedScalar::lookupOrAddToDict
(
"gamma2",
coeffDict_,
0.44
)
),
beta1_
(
dimensionedScalar::lookupOrAddToDict
(
"beta1",
coeffDict_,
0.075
)
),
beta2_
(
dimensionedScalar::lookupOrAddToDict
(
"beta2",
coeffDict_,
0.0828
)
),
betaStar_
(
dimensionedScalar::lookupOrAddToDict
(
"betaStar",
coeffDict_,
0.09
)
),
a1_
(
dimensionedScalar::lookupOrAddToDict
(
"a1",
coeffDict_,
0.31
)
),
b1_
(
dimensionedScalar::lookupOrAddToDict
(
"b1",
coeffDict_,
1.0
)
),
c1_
(
dimensionedScalar::lookupOrAddToDict
(
"c1",
coeffDict_,
10.0
)
),
F3_
(
Switch::lookupOrAddToDict
(
"F3",
coeffDict_,
false
)
),
y_(mesh_),
k_
(
IOobject
(
"k",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateIbK("k", mesh_, U_.db())
),
omega_
(
IOobject
(
"omega",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateIbOmega("omega", mesh_, U_.db())
),
nut_
(
IOobject
(
"nut",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateIbNut("nut", mesh_, U_.db())
)
{
bound(k_, k0_);
bound(omega_, omega0_);
nut_ =
(
a1_*k_/
max
(
a1_*omega_,
b1_*F23()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
)
);
nut_.correctBoundaryConditions();
printCoeffs();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<volSymmTensorField> immersedBoundaryKOmegaSST::R() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"R",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
k_.boundaryField().types()
)
);
}
tmp<volSymmTensorField> immersedBoundaryKOmegaSST::devReff() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"devRhoReff",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
-nuEff()*dev(twoSymm(fvc::grad(U_)))
)
);
}
tmp<fvVectorMatrix> immersedBoundaryKOmegaSST::divDevReff() const
{
const volScalarField nuEffective = nuEff();
return
(
- fvm::laplacian(nuEffective, U_)
- (fvc::grad(U_) & fvc::grad(nuEffective))
);
}
bool immersedBoundaryKOmegaSST::read()
{
if (RASModel::read())
{
alphaK1_.readIfPresent(coeffDict());
alphaK2_.readIfPresent(coeffDict());
alphaOmega1_.readIfPresent(coeffDict());
alphaOmega2_.readIfPresent(coeffDict());
gamma1_.readIfPresent(coeffDict());
gamma2_.readIfPresent(coeffDict());
beta1_.readIfPresent(coeffDict());
beta2_.readIfPresent(coeffDict());
betaStar_.readIfPresent(coeffDict());
a1_.readIfPresent(coeffDict());
b1_.readIfPresent(coeffDict());
c1_.readIfPresent(coeffDict());
F3_.readIfPresent("F3", coeffDict());
return true;
}
else
{
return false;
}
}
void immersedBoundaryKOmegaSST::correct()
{
// Bound in case of topological change
// HJ, 22/Aug/2007
if (mesh_.changing())
{
bound(k_, k0_);
bound(omega_, omega0_);
}
RASModel::correct();
if (!turbulence_)
{
return;
}
if (mesh_.changing())
{
y_.correct();
}
const volScalarField S2(2*magSqr(symm(fvc::grad(U_))));
volScalarField G("RASModel::G", nut_*S2);
// Update omega and G at the wall
omega_.boundaryField().updateCoeffs();
const volScalarField CDkOmega
(
(2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
);
const volScalarField F1(this->F1(CDkOmega));
// Turbulent frequency equation
fvScalarMatrix omegaEqn
(
fvm::ddt(omega_)
+ fvm::div(phi_, omega_)
+ fvm::SuSp(-fvc::div(phi_), omega_)
- fvm::laplacian(DomegaEff(F1), omega_)
==
gamma(F1)
*min
(
S2,
(c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2))
)
- fvm::Sp(beta(F1)*omega_, omega_)
- fvm::SuSp
(
(F1 - scalar(1))*CDkOmega/omega_,
omega_
)
);
omegaEqn.relax();
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// omegaEqn.completeAssembly();
solve(omegaEqn);
bound(omega_, omega0_);
// Turbulent kinetic energy equation
fvScalarMatrix kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
+ fvm::SuSp(-fvc::div(phi_), k_)
- fvm::laplacian(DkEff(F1), k_)
==
min(G, c1_*betaStar_*k_*omega_)
- fvm::Sp(betaStar_*omega_, k_)
);
kEqn.relax();
solve(kEqn);
bound(k_, k0_);
// Re-calculate viscosity
// Fixed sqrt(2) error. HJ, 10/Jun/2015
nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
nut_ = min(nut_, nuRatio()*nu());
nut_.correctBoundaryConditions();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,313 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / 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::incompressible::RASModels::immersedBoundaryKOmegaSST
Description
Implementation of the k-omega-SST turbulence model for incompressible
flows.
Turbulence model described in:
@verbatim
Menter, F., Esch, T.
"Elements of Industrial Heat Transfer Prediction"
16th Brazilian Congress of Mechanical Engineering (COBEM),
Nov. 2001.
@endverbatim
with updated coefficients from
@verbatim
Menter, F. R., Kuntz, M., and Langtry, R.,
"Ten Years of Industrial Experience with the SST Turbulence Model",
Turbulence, Heat and Mass Transfer 4, 2003,
pp. 625 - 632.
@endverbatim
but with the consistent production terms from the 2001 paper as form in the
2003 paper is a typo, see
@verbatim
http://turbmodels.larc.nasa.gov/sst.html
@endverbatim
and the addition of the optional F3 term for rough walls from
\verbatim
Hellsten, A.
"Some Improvements in Menters k-omega-SST turbulence model"
29th AIAA Fluid Dynamics Conference,
AIAA-98-2554,
June 1998.
\endverbatim
Note that this implementation is written in terms of alpha diffusion
coefficients rather than the more traditional sigma (alpha = 1/sigma) so
that the blending can be applied to all coefficuients in a consistent
manner. The paper suggests that sigma is blended but this would not be
consistent with the blending of the k-epsilon and k-omega models.
Also note that the error in the last term of equation (2) relating to
sigma has been corrected.
The default model coefficients correspond to the following:
@verbatim
immersedBoundaryKOmegaSSTCoeffs
{
alphaK1 0.85;
alphaK2 1.0;
alphaOmega1 0.5;
alphaOmega2 0.856;
beta1 0.075;
beta2 0.0828;
betaStar 0.09;
gamma1 5/9;
gamma2 0.44;
a1 0.31;
b1 1.0;
c1 10.0;
F3 no;
}
@endverbatim
SourceFiles
immersedBoundaryKOmegaSST.C
\*---------------------------------------------------------------------------*/
#ifndef immersedBoundaryKOmegaSST_H
#define immersedBoundaryKOmegaSST_H
#include "RASModel.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace RASModels
{
/*---------------------------------------------------------------------------*\
Class kOmega Declaration
\*---------------------------------------------------------------------------*/
class immersedBoundaryKOmegaSST
:
public RASModel
{
// Private data
// Model coefficients
dimensionedScalar alphaK1_;
dimensionedScalar alphaK2_;
dimensionedScalar alphaOmega1_;
dimensionedScalar alphaOmega2_;
dimensionedScalar gamma1_;
dimensionedScalar gamma2_;
dimensionedScalar beta1_;
dimensionedScalar beta2_;
dimensionedScalar betaStar_;
dimensionedScalar a1_;
dimensionedScalar b1_;
dimensionedScalar c1_;
Switch F3_;
//- Wall distance field
// Note: different to wall distance in parent RASModel
wallDist y_;
// Fields
volScalarField k_;
volScalarField omega_;
volScalarField nut_;
// Private member functions
tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
tmp<volScalarField> F2() const;
tmp<volScalarField> F3() const;
tmp<volScalarField> F23() const;
tmp<volScalarField> blend
(
const volScalarField& F1,
const dimensionedScalar& psi1,
const dimensionedScalar& psi2
) const
{
return F1*(psi1 - psi2) + psi2;
}
tmp<volScalarField> alphaK
(
const volScalarField& F1
) const
{
return blend(F1, alphaK1_, alphaK2_);
}
tmp<volScalarField> alphaOmega
(
const volScalarField& F1
) const
{
return blend(F1, alphaOmega1_, alphaOmega2_);
}
tmp<volScalarField> beta
(
const volScalarField& F1
) const
{
return blend(F1, beta1_, beta2_);
}
tmp<volScalarField> gamma
(
const volScalarField& F1
) const
{
return blend(F1, gamma1_, gamma2_);
}
public:
//- Runtime type information
TypeName("immersedBoundaryKOmegaSST");
// Constructors
//- Construct from components
immersedBoundaryKOmegaSST
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName = turbulenceModel::typeName,
const word& modelName = typeName
);
//- Destructor
virtual ~immersedBoundaryKOmegaSST()
{}
// Member Functions
//- Return the turbulence viscosity
virtual tmp<volScalarField> nut() const
{
return nut_;
}
//- Return the effective diffusivity for k
tmp<volScalarField> DkEff(const volScalarField& F1) const
{
return tmp<volScalarField>
(
new volScalarField("DkEff", alphaK(F1)*nut_ + nu())
);
}
//- Return the effective diffusivity for omega
tmp<volScalarField> DomegaEff(const volScalarField& F1) const
{
return tmp<volScalarField>
(
new volScalarField("DomegaEff", alphaOmega(F1)*nut_ + nu())
);
}
//- Return the turbulence kinetic energy
virtual tmp<volScalarField> k() const
{
return k_;
}
//- Return the turbulence specific dissipation rate
virtual tmp<volScalarField> omega() const
{
return omega_;
}
//- Return the turbulence kinetic energy dissipation rate
virtual tmp<volScalarField> epsilon() const
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"epsilon",
mesh_.time().timeName(),
mesh_
),
betaStar_*k_*omega_,
omega_.boundaryField().types()
)
);
}
//- Return the Reynolds stress tensor
virtual tmp<volSymmTensorField> R() const;
//- Return the effective stress tensor including the laminar stress
virtual tmp<volSymmTensorField> devReff() const;
//- Return the source term for the momentum equation
virtual tmp<fvVectorMatrix> divDevReff() const;
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct();
//- Read RASProperties dictionary
virtual bool read();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // namespace incompressible
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -91,6 +91,7 @@ immersedBoundaryKqRWallFunctionFvPatchField
)
{
this->checkType();
this->readPatchType(dict);
*this == this->patchInternalField();
}
@ -114,6 +115,7 @@ immersedBoundaryKqRWallFunctionFvPatchField
ptf.deadValue()
)
{
this->setPatchType(ptf);
this->checkType();
}
@ -133,6 +135,7 @@ immersedBoundaryKqRWallFunctionFvPatchField
ptf.deadValue()
)
{
this->setPatchType(ptf);
this->checkType();
}
@ -153,6 +156,7 @@ immersedBoundaryKqRWallFunctionFvPatchField
ptf.deadValue()
)
{
this->setPatchType(ptf);
this->checkType();
}

View file

@ -61,9 +61,11 @@ immersedBoundaryNutWallFunctionFvPatchScalarField
const dictionary& dict
)
:
nutkWallFunctionFvPatchScalarField(p, iF, dict),
nutkWallFunctionFvPatchScalarField(p, iF), // Do not read size
immersedBoundaryFieldBase<scalar>(p, true, 1e-6)
{}
{
this->readPatchType(dict);
}
immersedBoundaryNutWallFunctionFvPatchScalarField::
@ -77,7 +79,9 @@ immersedBoundaryNutWallFunctionFvPatchScalarField
:
nutkWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
immersedBoundaryFieldBase<scalar>(p, true, 1e-6)
{}
{
this->setPatchType(ptf);
}
immersedBoundaryNutWallFunctionFvPatchScalarField::
@ -88,7 +92,9 @@ immersedBoundaryNutWallFunctionFvPatchScalarField
:
nutkWallFunctionFvPatchScalarField(ptf),
immersedBoundaryFieldBase<scalar>(ptf.ibPatch(), true, 1e-6)
{}
{
this->setPatchType(ptf);
}
immersedBoundaryNutWallFunctionFvPatchScalarField::
@ -100,7 +106,9 @@ immersedBoundaryNutWallFunctionFvPatchScalarField
:
nutkWallFunctionFvPatchScalarField(ptf, iF),
immersedBoundaryFieldBase<scalar>(ptf.ibPatch(), true, 1e-6)
{}
{
this->setPatchType(ptf);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

View file

@ -60,14 +60,18 @@ immersedBoundaryOmegaWallFunctionFvPatchScalarField
const dictionary& dict
)
:
omegaWallFunctionFvPatchScalarField(p, iF, dict),
omegaWallFunctionFvPatchScalarField(p, iF), // Do not read size
immersedBoundaryFieldBase<scalar>
(
p,
Switch(dict.lookup("setDeadValue")),
readScalar(dict.lookup("deadValue"))
)
{}
{
this->readPatchType(dict);
*this == this->patchInternalField();
}
immersedBoundaryOmegaWallFunctionFvPatchScalarField::
@ -86,7 +90,9 @@ immersedBoundaryOmegaWallFunctionFvPatchScalarField
ptf.setDeadValue(),
ptf.deadValue()
)
{}
{
this->setPatchType(ptf);
}
immersedBoundaryOmegaWallFunctionFvPatchScalarField::
@ -102,7 +108,9 @@ immersedBoundaryOmegaWallFunctionFvPatchScalarField
ptf.setDeadValue(),
ptf.deadValue()
)
{}
{
this->setPatchType(ptf);
}
immersedBoundaryOmegaWallFunctionFvPatchScalarField::
@ -119,7 +127,9 @@ immersedBoundaryOmegaWallFunctionFvPatchScalarField
ptf.setDeadValue(),
ptf.deadValue()
)
{}
{
this->setPatchType(ptf);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -165,7 +175,7 @@ void immersedBoundaryOmegaWallFunctionFvPatchScalarField::updateCoeffs()
}
// Resize fields
if (size() != patch().size())
if (size() != this->ibPatch().size())
{
Info<< "Resizing immersedBoundaryOmegaWallFunction in evaluate: "
<< "patch: " << patch().size() << " field: " << size()
@ -191,7 +201,7 @@ void immersedBoundaryOmegaWallFunctionFvPatchScalarField::evaluate
)
{
// Resize fields
if (size() != patch().size())
if (size() != this->ibPatch().size())
{
Info<< "Resizing immersedBoundaryOmegaWallFunction in evaluate"
<< endl;