zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2009-01-20 13:35:56 +00:00
parent 21067168dd 5602a440d0
commit 372cd1b0e6
37 changed files with 1000 additions and 341 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
View File

@ -50,7 +50,7 @@ PDRkEpsilon::PDRkEpsilon
const volScalarField& rho,
const volVectorField& U,
const surfaceScalarField& phi,
basicThermo& thermophysicalModel
const basicThermo& thermophysicalModel
)
:
RASModel(typeName, rho, U, phi, thermophysicalModel),

2
applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H
View File

@ -93,7 +93,7 @@ public:
const volScalarField& rho,
const volVectorField& U,
const surfaceScalarField& phi,
basicThermo& thermophysicalModel
const basicThermo& thermophysicalModel
);

38
applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
View File

@ -129,9 +129,20 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
{
Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
<< "updateCoeffs() :"
<< " patch:" << patch().name()
<< " walltemperature "
<< " min:" << gMin(*this)
<< " max:" << gMax(*this)
<< " min:"
<< returnReduce
(
(this->size() > 0 ? min(*this) : VGREAT),
minOp<scalar>()
)
<< " max:"
<< returnReduce
(
(this->size() > 0 ? max(*this) : -VGREAT),
maxOp<scalar>()
)
<< " avg:" << gAverage(*this)
<< endl;
}
@ -163,7 +174,9 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
label nTotSize = returnReduce(this->size(), sumOp<label>());
Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
<< "updateCoeffs() : Out of " << nTotSize
<< "updateCoeffs() :"
<< " patch:" << patch().name()
<< " out of:" << nTotSize
<< " fixedBC:" << nFixed
<< " gradient:" << nTotSize-nFixed << endl;
}
@ -213,9 +226,22 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::evaluate
if (debug)
{
Info<< "Setting master and slave to wall temperature "
<< " min:" << gMin(*this)
<< " max:" << gMax(*this)
Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
<< "updateCoeffs() :"
<< " patch:" << patch().name()
<< " setting master and slave to wall temperature "
<< " min:"
<< returnReduce
(
(this->size() > 0 ? min(*this) : VGREAT),
minOp<scalar>()
)
<< " max:"
<< returnReduce
(
(this->size() > 0 ? max(*this) : -VGREAT),
maxOp<scalar>()
)
<< " avg:" << gAverage(*this)
<< endl;
}

2
applications/utilities/mesh/manipulation/cellSet/cellSetDict
View File

@ -82,7 +82,7 @@ topoSetSources
// Cells in cell zone
zoneToCell
{
name cellZone; // name of cellZone
name ".*Zone"; // name of cellZone, wildcards allowed.
}
// values of field within certain range

4
applications/utilities/mesh/manipulation/faceSet/faceSetDict
View File

@ -56,13 +56,13 @@ topoSetSources
// All faces of patch
patchToFace
{
name movingWall;
name ".*Wall"; // Name of patch, regular expressions allowed
}
// All faces of faceZone
zoneToFace
{
name faceZone1; // Name of faceZone
name ".*Zone1"; // Name of faceZone, regular expressions allowed
}
// Faces with face centre within box

6
applications/utilities/mesh/manipulation/pointSet/pointSetDict
View File

@ -57,6 +57,12 @@ topoSetSources
box (0 0 0) (1 1 1);
}
// All points in pointzone
zoneToPoint
{
name ".*Zone"; // name of pointZone, wildcards allowed.
}
// Select based on surface
surfaceToPoint
{

82
applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
View File

@ -488,17 +488,18 @@ labelList getNonRegionCells(const labelList& cellRegion, const label regionI)
}
// Get per region-region interface the sizes.
// If sumParallel does merge.
EdgeMap<label> getInterfaceSizes
// Get per region-region interface the sizes. If sumParallel sums sizes.
// Returns interfaces as straight list for looping in identical order.
void getInterfaceSizes
(
const polyMesh& mesh,
const labelList& cellRegion,
const bool sumParallel
const bool sumParallel,
edgeList& interfaces,
EdgeMap<label>& interfaceSizes
)
{
EdgeMap<label> interfaceSizes;
forAll(mesh.faceNeighbour(), faceI)
{
label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
@ -585,7 +586,12 @@ EdgeMap<label> getInterfaceSizes
}
}
return interfaceSizes;
// Make sure all processors have interfaces in same order
interfaces = interfaceSizes.toc();
if (sumParallel)
{
Pstream::scatter(interfaces);
}
}
@ -705,11 +711,7 @@ autoPtr<mapPolyMesh> createRegionMesh
if (otherRegion != -1)
{
edge interface
(
min(regionI, otherRegion),
max(regionI, otherRegion)
);
edge interface(regionI, otherRegion);
// Find the patch.
if (regionI < otherRegion)
@ -848,6 +850,7 @@ void createAndWriteRegion
const polyBoundaryMesh& newPatches = newMesh().boundaryMesh();
newPatches.checkParallelSync(true);
// Delete empty patches
// ~~~~~~~~~~~~~~~~~~~~
@ -863,13 +866,12 @@ void createAndWriteRegion
{
const polyPatch& pp = newPatches[patchI];
if
(
!isA<processorPolyPatch>(pp)
&& returnReduce(pp.size(), sumOp<label>()) > 0
)
if (!isA<processorPolyPatch>(pp))
{
oldToNew[patchI] = newI++;
if (returnReduce(pp.size(), sumOp<label>()) > 0)
{
oldToNew[patchI] = newI++;
}
}
}
@ -983,10 +985,15 @@ void createAndWriteRegion
}
// Create for every region-region interface with non-zero size two patches.
// First one is for minimumregion to maximumregion.
// Note that patches get created in same order on all processors (if parallel)
// since looping over synchronised 'interfaces'.
EdgeMap<label> addRegionPatches
(
fvMesh& mesh,
const regionSplit& cellRegion,
const edgeList& interfaces,
const EdgeMap<label>& interfaceSizes,
const wordList& regionNames
)
@ -998,15 +1005,12 @@ EdgeMap<label> addRegionPatches
EdgeMap<label> interfaceToPatch(cellRegion.nRegions());
// Keep start of added patches for later.
label minAddedPatchI = labelMax;
forAllConstIter(EdgeMap<label>, interfaceSizes, iter)
forAll(interfaces, interI)
{
if (iter() > 0)
{
const edge& e = iter.key();
const edge& e = interfaces[interI];
if (interfaceSizes[e] > 0)
{
label patchI = addPatch
(
mesh,
@ -1025,12 +1029,9 @@ EdgeMap<label> addRegionPatches
<< " " << mesh.boundaryMesh()[patchI].name()
<< endl;
interfaceToPatch.insert(iter.key(), patchI);
minAddedPatchI = min(minAddedPatchI, patchI);
interfaceToPatch.insert(e, patchI);
}
}
//Info<< "minAddedPatchI:" << minAddedPatchI << endl;
return interfaceToPatch;
}
@ -1348,24 +1349,26 @@ int main(int argc, char *argv[])
// Sizes of interface between regions. From pair of regions to number of
// faces.
EdgeMap<label> interfaceSizes
edgeList interfaces;
EdgeMap<label> interfaceSizes;
getInterfaceSizes
(
getInterfaceSizes
(
mesh,
cellRegion,
true // sum in parallel?
)
mesh,
cellRegion,
true, // sum in parallel?
interfaces,
interfaceSizes
);
Info<< "Region\tRegion\tFaces" << nl
<< "------\t------\t-----" << endl;
forAllConstIter(EdgeMap<label>, interfaceSizes, iter)
forAll(interfaces, interI)
{
const edge& e = iter.key();
const edge& e = interfaces[interI];
Info<< e[0] << '\t' << e[1] << '\t' << iter() << nl;
Info<< e[0] << '\t' << e[1] << '\t' << interfaceSizes[e] << nl;
}
Info<< endl;
@ -1511,6 +1514,7 @@ int main(int argc, char *argv[])
(
mesh,
cellRegion,
interfaces,
interfaceSizes,
regionNames
)

8
src/OSspecific/Unix/Unix.C
View File

@ -91,7 +91,9 @@ Foam::string Foam::getEnv(const word& envName)
}
else
{
return string::null;
// Return null-constructed string rather than string::null
// to avoid cyclic dependencies in the construction of globals
return string();
}
}
@ -277,7 +279,9 @@ Foam::fileName Foam::findEtcFile(const fileName& name, bool mandatory)
::exit(1);
}
return fileName::null;
// Return null-constructed fileName rather than fileName::null
// to avoid cyclic dependencies in the construction of globals
return fileName();
}

6
src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H
View File

@ -101,6 +101,12 @@ public:
//- Construct given size.
explicit inline DynamicList(const label);
//- Construct copy.
explicit inline DynamicList
(
const DynamicList<T, SizeInc, SizeMult, SizeDiv>&
);
//- Construct from UList. Size set to UList size.
explicit inline DynamicList(const UList<T>&);

11
src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H
View File

@ -49,6 +49,17 @@ inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList
}
template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList
(
const DynamicList<T, SizeInc, SizeMult, SizeDiv>& lst
)
:
List<T>(lst),
capacity_(lst.size())
{}
template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList
(

5
src/OpenFOAM/global/global.Cver
View File

@ -38,11 +38,6 @@ Description
const char* const Foam::FOAMversion = "VERSION_STRING";
const char* const Foam::FOAMbuild = "BUILD_STRING";
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Static initializers for string::null, word::null and fileName::null
#include "stringsGlobals.C"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Setup an error handler for the global new operator

4
src/OpenFOAM/include/OSspecific.H
View File

@ -62,7 +62,7 @@ pid_t pgid();
bool env(const word&);
//- Return environment variable of given name
// Return string::null if the environment is undefined
// Return string() if the environment is undefined
string getEnv(const word& name);
//- Set an environment variable
@ -101,7 +101,7 @@ bool chDir(const fileName& dir);
// -# shipped settings:
// - $WM_PROJECT_DIR/etc/
//
// @return the full path name or fileName::null if the name cannot be found
// @return the full path name or fileName() if the name cannot be found
// Optionally abort if the file cannot be found
fileName findEtcFile(const fileName& name, bool mandatory=false);

1
src/OpenFOAM/primitives/strings/fileName/fileName.C
View File

@ -33,6 +33,7 @@ License
const char* const Foam::fileName::typeName = "fileName";
int Foam::fileName::debug(debug::debugSwitch(fileName::typeName, 0));
const Foam::fileName Foam::fileName::null;
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

3
src/OpenFOAM/primitives/strings/string/string.C
View File

@ -29,6 +29,9 @@ License
/* * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * */
const char* const Foam::string::typeName = "string";
int Foam::string::debug(debug::debugSwitch(string::typeName, 0));
const Foam::string Foam::string::null;
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

8
src/OpenFOAM/primitives/strings/string/string.H
View File

@ -80,6 +80,8 @@ public:
// Static data members
static const char* const typeName;
static int debug;
static const string null;
@ -116,8 +118,6 @@ public:
// Member Functions
// Access
//- Count and return the number of a given character in the string
size_type count(const char) const;
@ -130,9 +130,6 @@ public:
template<class String>
static inline bool meta(const string&, const char quote='\\');
// Edit
//- Strip invalid characters from the given string
template<class String>
static inline bool stripInvalid(string&);
@ -145,7 +142,6 @@ public:
template<class String>
static inline string quotemeta(const string&, const char quote='\\');
//- Replace first occurence of sub-string oldStr with newStr
// starting at start
string& replace

4
src/OpenFOAM/primitives/strings/stringsGlobals.C
View File

@ -38,8 +38,4 @@ Description
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::string Foam::string::null;
const Foam::word Foam::word::null;
const Foam::fileName Foam::fileName::null;
// ************************************************************************* //

1
src/OpenFOAM/primitives/strings/word/word.C
View File

@ -31,5 +31,6 @@ License
const char* const Foam::word::typeName = "word";
int Foam::word::debug(Foam::debug::debugSwitch(word::typeName, 0));
const Foam::word Foam::word::null;
// ************************************************************************* //

15
src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
View File

@ -211,7 +211,6 @@ Foam::refinementSurfaces::refinementSurfaces
minLevel_[globalRegionI] = iter();
maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
// Check validity
if
@ -231,6 +230,13 @@ Foam::refinementSurfaces::refinementSurfaces
<< exit(FatalError);
}
}
forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
{
label globalRegionI = regionOffset_[surfI] + iter.key();
perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
}
//// Optional patch names and patch types
//forAllConstIter(HashTable<word>, regionPatchName[surfI], iter)
@ -387,7 +393,6 @@ Foam::refinementSurfaces::refinementSurfaces
minLevel_[globalRegionI] = iter();
maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
// Check validity
if
@ -407,6 +412,12 @@ Foam::refinementSurfaces::refinementSurfaces
<< exit(FatalError);
}
}
forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
{
label globalRegionI = regionOffset_[surfI] + iter.key();
perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
}
}
}

1
src/meshTools/Make/files
View File

@ -56,6 +56,7 @@ indexedOctree/treeDataTriSurface.C
searchableSurface = searchableSurface
$(searchableSurface)/distributedTriSurfaceMesh.C
$(searchableSurface)/searchableBox.C
$(searchableSurface)/searchablePlane.C
$(searchableSurface)/searchablePlate.C
$(searchableSurface)/searchableSphere.C
$(searchableSurface)/searchableSurface.C

231
src/meshTools/searchableSurface/searchablePlane.C Normal file
View File

@ -0,0 +1,231 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "searchablePlane.H"
#include "addToRunTimeSelectionTable.H"
#include "SortableList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(searchablePlane, 0);
addToRunTimeSelectionTable(searchableSurface, searchablePlane, dict);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::pointIndexHit Foam::searchablePlane::findLine
(
const point& start,
const point& end
) const
{
pointIndexHit info(true, vector::zero, 0);
linePointRef l(start, end);
scalar t = lineIntersect(l);
if (t < 0 || t > 1)
{
info.setMiss();
info.setIndex(-1);
}
else
{
info.setPoint(start+t*l.vec());
}
return info;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchablePlane::searchablePlane
(
const IOobject& io,
const point& basePoint,
const vector& normal
)
:
searchableSurface(io),
plane(basePoint, normal)
{}
Foam::searchablePlane::searchablePlane
(
const IOobject& io,
const dictionary& dict
)
:
searchableSurface(io),
plane(dict)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::searchablePlane::~searchablePlane()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::wordList& Foam::searchablePlane::regions() const
{
if (regions_.empty())
{
regions_.setSize(1);
regions_[0] = "region0";
}
return regions_;
}
void Foam::searchablePlane::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
List<pointIndexHit>& info
) const
{
info.setSize(samples.size());
forAll(samples, i)
{
info[i].setPoint(nearestPoint(samples[i]));
if (magSqr(samples[i]-info[i].rawPoint()) > nearestDistSqr[i])
{
info[i].setIndex(-1);
info[i].setMiss();
}
else
{
info[i].setIndex(0);
info[i].setHit();
}
}
}
void Foam::searchablePlane::findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
info.setSize(start.size());
forAll(start, i)
{
info[i] = findLine(start[i], end[i]);
}
}
void Foam::searchablePlane::findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
findLine(start, end, info);
}
void Foam::searchablePlane::findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >& info
) const
{
List<pointIndexHit> nearestInfo;
findLine(start, end, nearestInfo);
info.setSize(start.size());
forAll(info, pointI)
{
if (nearestInfo[pointI].hit())
{
info[pointI].setSize(1);
info[pointI][0] = nearestInfo[pointI];
}
else
{
info[pointI].clear();
}
}
}
void Foam::searchablePlane::getRegion
(
const List<pointIndexHit>& info,
labelList& region
) const
{
region.setSize(info.size());
region = 0;
}
void Foam::searchablePlane::getNormal
(
const List<pointIndexHit>& info,
vectorField& n
) const
{
n.setSize(info.size());
n = normal();
}
void Foam::searchablePlane::getVolumeType
(
const pointField& points,
List<volumeType>& volType
) const
{
FatalErrorIn
(
"searchableCollection::getVolumeType(const pointField&"
", List<volumeType>&) const"
) << "Volume type not supported for plane."
<< exit(FatalError);
}
// ************************************************************************* //

199
src/meshTools/searchableSurface/searchablePlane.H Normal file
View File

@ -0,0 +1,199 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::searchablePlane
Description
Searching on plane. See plane.H
SourceFiles
searchablePlane.C
\*---------------------------------------------------------------------------*/
#ifndef searchablePlane_H
#define searchablePlane_H
#include "searchableSurface.H"
#include "plane.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
/*---------------------------------------------------------------------------*\
Class searchablePlane Declaration
\*---------------------------------------------------------------------------*/
class searchablePlane
:
public searchableSurface,
public plane
{
private:
// Private Member Data
mutable wordList regions_;
// Private Member Functions
pointIndexHit findLine
(
const point& start,
const point& end
) const;
//- Disallow default bitwise copy construct
searchablePlane(const searchablePlane&);
//- Disallow default bitwise assignment
void operator=(const searchablePlane&);
public:
//- Runtime type information
TypeName("searchablePlane");
// Constructors
//- Construct from components
searchablePlane
(
const IOobject& io,
const point& basePoint,
const vector& normal
);
//- Construct from dictionary (used by searchableSurface)
searchablePlane
(
const IOobject& io,
const dictionary& dict
);
// Destructor
virtual ~searchablePlane();
// Member Functions
virtual const wordList& regions() const;
//- Whether supports volume type below
virtual bool hasVolumeType() const
{
return false;
}
//- Range of local indices that can be returned.
virtual label size() const
{
return 1;
}
// Multiple point queries.
virtual void findNearest
(
const pointField& sample,
const scalarField& nearestDistSqr,
List<pointIndexHit>&
) const;
virtual void findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
virtual void findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
//- Get all intersections in order from start to end.
virtual void findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >&
) const;
//- From a set of points and indices get the region
virtual void getRegion
(
const List<pointIndexHit>&,
labelList& region
) const;
//- From a set of points and indices get the normal
virtual void getNormal
(
const List<pointIndexHit>&,
vectorField& normal
) const;
//- Determine type (inside/outside/mixed) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
virtual void getVolumeType
(
const pointField&,
List<volumeType>&
) const;
// regIOobject implementation
bool writeData(Ostream&) const
{
notImplemented("searchablePlane::writeData(Ostream&) const");
return false;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

22
src/meshTools/searchableSurface/searchablePlate.C
View File

@ -218,7 +218,16 @@ Foam::searchablePlate::searchablePlate
origin_(origin),
span_(span),
normalDir_(calcNormal(span_))
{}
{
if (debug)
{
Info<< "searchablePlate::searchablePlate :"
<< " origin:" << origin_
<< " origin+span:" << origin_+span_
<< " normal:" << vector::componentNames[normalDir_]
<< endl;
}
}
Foam::searchablePlate::searchablePlate
@ -231,7 +240,16 @@ Foam::searchablePlate::searchablePlate
origin_(dict.lookup("origin")),
span_(dict.lookup("span")),
normalDir_(calcNormal(span_))
{}
{
if (debug)
{
Info<< "searchablePlate::searchablePlate :"
<< " origin:" << origin_
<< " origin+span:" << origin_+span_
<< " normal:" << vector::componentNames[normalDir_]
<< endl;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //

9
src/meshTools/searchableSurface/searchablePlate.H
View File

@ -27,7 +27,14 @@ Class
Description
Searching on finite plate. Plate has to be aligned with coordinate
axes!
axes.
Plate defined as origin and span. One of the components of span has
to be 0 which defines the normal direction. E.g.
span = (Sx Sy 0) // plate in x-y plane
origin = (Ox Oy Oz)
now plane is from (Ox Oy Oz) to (Ox+Sx Oy+Sy Oz)
SourceFiles
searchablePlate.C

38
src/meshTools/sets/cellSources/nbrToCell/nbrToCell.C
View File

@ -22,8 +22,6 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
\*---------------------------------------------------------------------------*/
#include "nbrToCell.H"
@ -58,22 +56,46 @@ Foam::topoSetSource::addToUsageTable Foam::nbrToCell::usage_
void Foam::nbrToCell::combine(topoSet& set, const bool add) const
{
const cellList& cells = mesh().cells();
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
boolList isCoupled(mesh_.nFaces()-mesh_.nInternalFaces(), false);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
label faceI = pp.start();
forAll(pp, i)
{
isCoupled[faceI-mesh_.nInternalFaces()] = true;
faceI++;
}
}
}
forAll(cells, cellI)
{
const cell& cll = cells[cellI];
const cell& cFaces = cells[cellI];
label nInternalFaces = 0;
label nNbrCells = 0;
forAll(cll, i)
forAll(cFaces, i)
{
if (mesh().isInternalFace(cll[i]))
label faceI = cFaces[i];
if (mesh_.isInternalFace(faceI))
{
nInternalFaces++;
nNbrCells++;
}
else if (isCoupled[faceI-mesh_.nInternalFaces()])
{
nNbrCells++;
}
}
if (nInternalFaces <= minNbrs_)
if (nNbrCells <= minNbrs_)
{
addOrDelete(set, cellI, add);
}

2
src/meshTools/sets/cellSources/nbrToCell/nbrToCell.H
View File

@ -27,7 +27,7 @@ Class
Description
A topoSetSource to select cells based on number of neighbouring cells
(i.e. number of internal faces)
(i.e. number of internal or coupled faces)
SourceFiles
nbrToCell.C

81
src/meshTools/sets/faceSources/cellToFace/cellToFace.C
View File

@ -22,14 +22,13 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
\*---------------------------------------------------------------------------*/
#include "cellToFace.H"
#include "polyMesh.H"
#include "cellSet.H"
#include "Time.H"
#include "syncTools.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -102,48 +101,58 @@ void Foam::cellToFace::combine(topoSet& set, const bool add) const
{
// Add all faces whose both neighbours are in set.
// Count number of cells using face.
Map<label> numCells(loadedSet.size());
label nInt = mesh_.nInternalFaces();
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
for
(
cellSet::const_iterator iter = loadedSet.begin();
iter != loadedSet.end();
++iter
)
// Check all internal faces
for (label faceI = 0; faceI < nInt; faceI++)
{
label cellI = iter.key();
const labelList& cFaces = mesh_.cells()[cellI];
forAll(cFaces, cFaceI)
if (loadedSet.found(own[faceI]) && loadedSet.found(nei[faceI]))
{
label faceI = cFaces[cFaceI];
Map<label>::iterator fndFace = numCells.find(faceI);
if (fndFace == numCells.end())
{
numCells.insert(faceI, 1);
}
else
{
fndFace()++;
}
addOrDelete(set, faceI, add);
}
}
// Include faces that are referenced twice
for
(
Map<label>::const_iterator iter = numCells.begin();
iter != numCells.end();
++iter
)
// Get coupled cell status
boolList neiInSet(mesh_.nFaces()-nInt, false);
forAll(patches, patchI)
{
if (iter() == 2)
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
addOrDelete(set, iter.key(), add);
label faceI = pp.start();
forAll(pp, i)
{
neiInSet[faceI-nInt] = loadedSet.found(own[faceI]);
faceI++;
}
}
}
syncTools::swapBoundaryFaceList(mesh_, neiInSet, false);
// Check all boundary faces
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
label faceI = pp.start();
forAll(pp, i)
{
if (loadedSet.found(own[faceI]) && neiInSet[faceI-nInt])
{
addOrDelete(set, faceI, add);
}
faceI++;
}
}
}
}

1
src/meshTools/sets/faceSources/cellToFace/cellToFace.H
View File

@ -30,6 +30,7 @@ Description
Either all faces of cell or some other criterion.
See implementation.
Note: when picking up coupled faces uses cells on neighbouring processors.
SourceFiles
cellToFace.C

86
src/sampling/sampledSurface/distanceSurface/distanceSurface.C
View File

@ -45,6 +45,11 @@ namespace Foam
void Foam::distanceSurface::createGeometry()
{
if (debug)
{
Pout<< "distanceSurface::createGeometry :updating geometry." << endl;
}
// Clear any stored topologies
facesPtr_.clear();
@ -67,7 +72,7 @@ void Foam::distanceSurface::createGeometry()
false
),
fvm,
dimensionedScalar("zero", dimless/dimTime, 0)
dimensionedScalar("zero", dimLength, 0)
)
);
volScalarField& cellDistance = cellDistancePtr_();
@ -157,6 +162,7 @@ void Foam::distanceSurface::createGeometry()
}
}
// On processor patches the mesh.C() will already be the cell centre
// on the opposite side so no need to swap cellDistance.
@ -164,23 +170,70 @@ void Foam::distanceSurface::createGeometry()
// Distance to points
pointDistance_.setSize(fvm.nPoints());
{
const pointField& pts = fvm.points();
List<pointIndexHit> nearest;
surfPtr_().findNearest
(
fvm.points(),
scalarField(fvm.nPoints(), GREAT),
pts,
scalarField(pts.size(), GREAT),
nearest
);
forAll(pointDistance_, pointI)
if (signed_)
{
pointDistance_[pointI] = Foam::mag
(
nearest[pointI].hitPoint()
- fvm.points()[pointI]
);
vectorField normal;
surfPtr_().getNormal(nearest, normal);
forAll(nearest, i)
{
vector d(pts[i]-nearest[i].hitPoint());
if ((d&normal[i]) > 0)
{
pointDistance_[i] = Foam::mag(d);
}
else
{
pointDistance_[i] = -Foam::mag(d);
}
}
}
else
{
forAll(nearest, i)
{
pointDistance_[i] = Foam::mag(pts[i]-nearest[i].hitPoint());
}
}
}
if (debug)
{
Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
cellDistance.write();
pointScalarField pDist
(
IOobject
(
"pointDistance",
fvm.time().timeName(),
fvm.time(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
pointMesh::New(fvm),
dimensionedScalar("zero", dimLength, 0)
);
pDist.internalField() = pointDistance_;
Pout<< "Writing point distance:" << pDist.objectPath() << endl;
pDist.write();
}
//- Direct from cell field and point field.
isoSurfPtr_.reset
(
@ -196,6 +249,7 @@ void Foam::distanceSurface::createGeometry()
if (debug)
{
print(Pout);
Pout<< endl;
}
}
@ -264,6 +318,13 @@ bool Foam::distanceSurface::needsUpdate() const
bool Foam::distanceSurface::expire()
{
if (debug)
{
Pout<< "distanceSurface::expire :"
<< " have-facesPtr_:" << facesPtr_.valid()
<< " needsUpdate_:" << needsUpdate_ << endl;
}
// Clear any stored topologies
facesPtr_.clear();
@ -280,6 +341,13 @@ bool Foam::distanceSurface::expire()
bool Foam::distanceSurface::update()
{
if (debug)
{
Pout<< "distanceSurface::update :"
<< " have-facesPtr_:" << facesPtr_.valid()
<< " needsUpdate_:" << needsUpdate_ << endl;
}
if (!needsUpdate_)
{
return false;

97
src/sampling/sampledSurface/isoSurface/isoSurface.C
View File

@ -58,6 +58,31 @@ Foam::scalar Foam::isoSurface::isoFraction
}
bool Foam::isoSurface::isEdgeOfFaceCut
(
const scalarField& pVals,
const face& f,
const bool ownLower,
const bool neiLower
) const
{
forAll(f, fp)
{
bool fpLower = (pVals[f[fp]] < iso_);
if
(
(fpLower != ownLower)
|| (fpLower != neiLower)
|| (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
)
{
return true;
}
}
return false;
}
// Determine for every face/cell whether it (possibly) generates triangles.
void Foam::isoSurface::calcCutTypes
(
@ -84,22 +109,13 @@ void Foam::isoSurface::calcCutTypes
}
else
{
// Mesh edge.
// See if any mesh edge is cut by looping over all the edges of the
// face.
const face f = mesh_.faces()[faceI];
forAll(f, fp)
if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
{
bool fpLower = (pVals[f[fp]] < iso_);
if
(
(fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
|| (fpLower != ownLower)
|| (fpLower != neiLower)
)
{
faceCutType_[faceI] = CUT;
break;
}
faceCutType_[faceI] = CUT;
}
}
}
@ -117,22 +133,13 @@ void Foam::isoSurface::calcCutTypes
{
bool ownLower = (cVals[own[faceI]] < iso_);
// Mesh edge.
const face f = mesh_.faces()[faceI];
forAll(f, fp)
if (isEdgeOfFaceCut(pVals, f, ownLower, ownLower))
{
bool fpLower = (pVals[f[fp]] < iso_);
if
(
(fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
|| (fpLower != ownLower)
)
{
faceCutType_[faceI] = CUT;
break;
}
faceCutType_[faceI] = CUT;
}
faceI++;
}
}
@ -152,19 +159,9 @@ void Foam::isoSurface::calcCutTypes
// Mesh edge.
const face f = mesh_.faces()[faceI];
forAll(f, fp)
if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
{
bool fpLower = (pVals[f[fp]] < iso_);
if
(
(fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
|| (fpLower != ownLower)
|| (fpLower != neiLower)
)
{
faceCutType_[faceI] = CUT;
break;
}
faceCutType_[faceI] = CUT;
}
}
faceI++;
@ -1355,6 +1352,30 @@ Foam::isoSurface::isoSurface
iso_(iso),
mergeDistance_(mergeTol*mesh_.bounds().mag())
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
const labelList& own = mesh_.faceOwner();
// Check
forAll(patches, patchI)
{
if (isA<emptyPolyPatch>(patches[patchI]))
{
FatalErrorIn
(
"isoSurface::isoSurface\n"
"(\n"
" const volScalarField& cVals,\n"
" const scalarField& pVals,\n"
" const scalar iso,\n"
" const bool regularise,\n"
" const scalar mergeTol\n"
")\n"
) << "Iso surfaces not supported on case with empty patches."
<< exit(FatalError);
}
}
// Determine if any cut through face/cell
calcCutTypes(cVals, pVals);
@ -1364,8 +1385,6 @@ Foam::isoSurface::isoSurface
PackedList<1> isBoundaryPoint(mesh_.nPoints());
labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
const labelList& own = mesh_.faceOwner();
forAll(patches, patchI)
{

14
src/sampling/sampledSurface/isoSurface/isoSurface.H
View File

@ -31,6 +31,9 @@ Description
G.M. Treece, R.W. Prager and A.H. Gee.
Note:
- not possible on patches of type 'empty'. There are no values on
'empty' patch fields so even the api would have to change
(no volScalarField as argument). Too messy.
- in parallel the regularisation (coarsening) always takes place
and slightly different surfaces will be created compared to non-parallel.
The surface will still be continuous though!
@ -123,6 +126,15 @@ class isoSurface
const scalar s1
) const;
//- Check if any edge of a face is cut
bool isEdgeOfFaceCut
(
const scalarField& pVals,
const face& f,
const bool ownLower,
const bool neiLower
) const;
//- Set faceCutType,cellCutType.
void calcCutTypes
(
@ -217,7 +229,7 @@ class isoSurface
) const;
template<class Type>
label generateTriPoints
label generateFaceTriPoints
(
const volScalarField& cVals,
const scalarField& pVals,

25
src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
View File

@ -200,7 +200,7 @@ void Foam::isoSurface::generateTriPoints
template<class Type>
Foam::label Foam::isoSurface::generateTriPoints
Foam::label Foam::isoSurface::generateFaceTriPoints
(
const volScalarField& cVals,
const scalarField& pVals,
@ -288,6 +288,19 @@ void Foam::isoSurface::generateTriPoints
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
if
(
(cVals.size() != mesh_.nCells())
|| (pVals.size() != mesh_.nPoints())
|| (cCoords.size() != mesh_.nCells())
|| (pCoords.size() != mesh_.nPoints())
|| (snappedCc.size() != mesh_.nCells())
|| (snappedPoint.size() != mesh_.nPoints())
)
{
FatalErrorIn("isoSurface::generateTriPoints(..)")
<< "Incorrect size." << abort(FatalError);
}
// Determine neighbouring snap status
labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
@ -319,7 +332,7 @@ void Foam::isoSurface::generateTriPoints
{
if (faceCutType_[faceI] != NOTCUT)
{
generateTriPoints
generateFaceTriPoints
(
cVals,
pVals,
@ -357,7 +370,7 @@ void Foam::isoSurface::generateTriPoints
{
if (faceCutType_[faceI] != NOTCUT)
{
generateTriPoints
generateFaceTriPoints
(
cVals,
pVals,
@ -384,14 +397,14 @@ void Foam::isoSurface::generateTriPoints
}
else if (isA<emptyPolyPatch>(pp))
{
// Assume zero-gradient.
// Assume zero-gradient. But what about coordinates?
label faceI = pp.start();
forAll(pp, i)
{
if (faceCutType_[faceI] != NOTCUT)
{
generateTriPoints
generateFaceTriPoints
(
cVals,
pVals,
@ -423,7 +436,7 @@ void Foam::isoSurface::generateTriPoints
{
if (faceCutType_[faceI] != NOTCUT)
{
generateTriPoints
generateFaceTriPoints
(
cVals,
pVals,

103
src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
View File

@ -26,7 +26,6 @@ License
#include "SpalartAllmaras.H"
#include "addToRunTimeSelectionTable.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -54,7 +53,7 @@ tmp<volScalarField> SpalartAllmaras::fv1() const
tmp<volScalarField> SpalartAllmaras::fv2() const
{
volScalarField chi = nuTilda_/nu();
return 1.0/pow3(scalar(1) + chi/Cv2_);
return 1/pow3(scalar(1) + chi/Cv2_);
}
@ -71,70 +70,68 @@ tmp<volScalarField> SpalartAllmaras::fv3() const
}
tmp<volScalarField> SpalartAllmaras::calcS(const volTensorField& gradU)
tmp<volScalarField> SpalartAllmaras::S(const volTensorField& gradU) const
{
return ::sqrt(2.0)*mag(skew(gradU));
return sqrt(2.0)*mag(skew(gradU));
}
tmp<volScalarField> SpalartAllmaras::calcSTilda(const volTensorField& gradU)
tmp<volScalarField> SpalartAllmaras::STilda
(
const volScalarField& S,
const volScalarField& dTilda
) const
{
return fv3()*calcS(gradU) + fv2()*nuTilda_/sqr(kappa_*dTilda_);
return fv3()*S + fv2()*nuTilda_/sqr(kappa_*dTilda);
}
tmp<volScalarField> SpalartAllmaras::r
(
const volScalarField& visc,
const volScalarField& S
const volScalarField& S,
const volScalarField& dTilda
) const
{
tmp<volScalarField> tr
return min
(
new volScalarField
(
min
(
visc
/(
max
(
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)
*sqr(kappa_*dTilda_)
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
),
scalar(10.0)
)
)
visc
/(
max
(
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)
*sqr(kappa_*dTilda)
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
),
scalar(10)
);
return tr;
}
tmp<volScalarField> SpalartAllmaras::fw(const volScalarField& S) const
tmp<volScalarField> SpalartAllmaras::fw
(
const volScalarField& S,
const volScalarField& dTilda
) const
{
volScalarField r = this->r(nuTilda_, S);
volScalarField r = this->r(nuTilda_, S, dTilda);
volScalarField g = r + Cw2_*(pow6(r) - r);
return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
}
void SpalartAllmaras::dTildaUpdate(const volScalarField&)
tmp<volScalarField> SpalartAllmaras::dTilda(const volScalarField&) const
{
if (mesh_.changing())
{
dTilda_ = min(CDES_*delta(), wallDist(mesh_).y());
}
return min(CDES_*delta(), y_);
}
@ -242,6 +239,8 @@ SpalartAllmaras::SpalartAllmaras
)
),
y_(mesh_),
nuTilda_
(
IOobject
@ -255,8 +254,6 @@ SpalartAllmaras::SpalartAllmaras
mesh_
),
dTilda_(min(CDES_*delta(), wallDist(mesh_).y())),
nuSgs_
(
IOobject
@ -277,11 +274,15 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
{
LESModel::correct(gradU);
const volScalarField STilda = calcSTilda(gradU);
if (mesh_.changing())
{
y_.correct();
y_.boundaryField() = max(y_.boundaryField(), VSMALL);
}
const volScalarField S = calcS(gradU);
dTildaUpdate(S);
const volScalarField S = this->S(gradU);
const volScalarField dTilda = this->dTilda(S);
const volScalarField STilda = this->STilda(S, dTilda);
fvScalarMatrix nuTildaEqn
(
@ -296,7 +297,7 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
- alphaNut_*Cb2_*magSqr(fvc::grad(nuTilda_))
==
Cb1_*STilda*nuTilda_
- fvm::Sp(Cw1_*fw(STilda)*nuTilda_/sqr(dTilda_), nuTilda_)
- fvm::Sp(Cw1_*fw(STilda, dTilda)*nuTilda_/sqr(dTilda), nuTilda_)
);
nuTildaEqn.relax();
@ -310,9 +311,15 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
}
tmp<volScalarField> SpalartAllmaras::k() const
{
return sqr(nuSgs()/ck_/dTilda(S(fvc::grad(U()))));
}
tmp<volScalarField> SpalartAllmaras::epsilon() const
{
return 2.0*nuEff()*magSqr(symm(fvc::grad(U())));
return 2*nuEff()*magSqr(symm(fvc::grad(U())));
}

47
src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
View File

@ -38,6 +38,7 @@ SourceFiles
#include "LESModel.H"
#include "volFields.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -86,30 +87,39 @@ protected:
// Fields
wallDist y_;
volScalarField nuTilda_;
volScalarField dTilda_;
volScalarField nuSgs_;
// Protected member functions
// Helper functions
virtual tmp<volScalarField> fv1() const;
virtual tmp<volScalarField> fv2() const;
virtual tmp<volScalarField> fv3() const;
virtual tmp<volScalarField> S(const volTensorField& gradU) const;
virtual tmp<volScalarField> fv1() const;
virtual tmp<volScalarField> fv2() const;
virtual tmp<volScalarField> fv3() const;
//-
virtual tmp<volScalarField> calcS(const volTensorField& gradU);
virtual tmp<volScalarField> calcSTilda(const volTensorField& gradU);
virtual tmp<volScalarField> r
(
const volScalarField& visc,
const volScalarField& S
) const;
virtual tmp<volScalarField> fw(const volScalarField& S) const;
virtual tmp<volScalarField> STilda
(
const volScalarField& S,
const volScalarField& dTilda
) const;
//- Length scale calculation
virtual void dTildaUpdate(const volScalarField& S);
virtual tmp<volScalarField> r
(
const volScalarField& visc,
const volScalarField& S,
const volScalarField& dTilda
) const;
virtual tmp<volScalarField> fw
(
const volScalarField& S,
const volScalarField& dTilda
) const;
//- Length scale
virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
public:
@ -138,10 +148,7 @@ public:
// Member Functions
//- Return SGS kinetic energy
virtual tmp<volScalarField> k() const
{
return sqr(nuSgs()/ck_/dTilda_);
}
virtual tmp<volScalarField> k() const;
//- Return sub-grid disipation rate
virtual tmp<volScalarField> epsilon() const;

58
src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.C
View File

@ -26,7 +26,6 @@ License
#include "SpalartAllmarasDDES.H"
#include "addToRunTimeSelectionTable.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -50,51 +49,42 @@ tmp<volScalarField> SpalartAllmarasDDES::rd
const volScalarField& S
) const
{
volScalarField d = wallDist(mesh_).y();
tmp<volScalarField> trd
return min
(
new volScalarField
(
min
visc
/(
max
(
visc
/(
max
(
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)*sqr(kappa_*d)
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
), scalar(10.0)
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)*sqr(kappa_*y_)
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
)
),
scalar(10)
);
return trd;
}
tmp<volScalarField> SpalartAllmarasDDES::fd(const volScalarField& S)
tmp<volScalarField> SpalartAllmarasDDES::fd(const volScalarField& S) const
{
return 1.0 - tanh(pow3(8.0*rd(nuSgs_ + transport().nu(), S)));
return 1 - tanh(pow3(8*rd(nuEff(), S)));
}
void SpalartAllmarasDDES::dTildaUpdate(const volScalarField& S)
tmp<volScalarField> SpalartAllmarasDDES::dTilda(const volScalarField& S) const
{
dTilda_ =
wallDist(mesh_).y()
- fd(S)*max
(
dimensionedScalar("zero", dimLength, 0.0),
wallDist(mesh_).y() - CDES_*delta()
);
return max
(
y_
- fd(S)
*max(y_ - CDES_*delta(), dimensionedScalar("zero", dimLength, 0)),
dimensionedScalar("small", dimLength, SMALL)
);
}

19
src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.H
View File

@ -62,6 +62,14 @@ class SpalartAllmarasDDES
{
// Private member functions
tmp<volScalarField> fd(const volScalarField& S) const;
tmp<volScalarField> rd
(
const volScalarField& visc,
const volScalarField& S
) const;
// Disallow default bitwise copy construct and assignment
SpalartAllmarasDDES(const SpalartAllmarasDDES&);
SpalartAllmarasDDES& operator=(const SpalartAllmarasDDES&);
@ -71,15 +79,8 @@ protected:
// Protected member functions
tmp<volScalarField> fd(const volScalarField& S);
tmp<volScalarField> rd
(
const volScalarField& visc,
const volScalarField& S
) const;
//- Length scale calculation
virtual void dTildaUpdate(const volScalarField& S);
//- Length scale
virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
public:

96
src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C
View File

@ -26,7 +26,6 @@ License
#include "SpalartAllmarasIDDES.H"
#include "addToRunTimeSelectionTable.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -46,22 +45,29 @@ addToRunTimeSelectionTable(LESModel, SpalartAllmarasIDDES, dictionary);
tmp<volScalarField> SpalartAllmarasIDDES::alpha() const
{
return
0.25
- wallDist(mesh_).y()
/dimensionedScalar("hMax", dimLength, max(cmptMax(delta())));
return max
(
0.25 - y_/dimensionedScalar("hMax", dimLength, max(cmptMax(delta()))),
scalar(-5)
);
}
tmp<volScalarField> SpalartAllmarasIDDES::ft(const volScalarField& S) const
tmp<volScalarField> SpalartAllmarasIDDES::ft
(
const volScalarField& S
) const
{
return tanh(pow3(sqr(ct_)*r(nuSgs_, S)));
return tanh(pow3(sqr(ct_)*rd(nuSgs_, S)));
}
tmp<volScalarField> SpalartAllmarasIDDES::fl(const volScalarField& S) const
tmp<volScalarField> SpalartAllmarasIDDES::fl
(
const volScalarField& S
) const
{
return tanh(pow(sqr(cl_)*r(transport().nu(), S), 10));
return tanh(pow(sqr(cl_)*rd(nu(), S), 10));
}
@ -71,33 +77,24 @@ tmp<volScalarField> SpalartAllmarasIDDES::rd
const volScalarField& S
) const
{
volScalarField d = wallDist(mesh_).y();
tmp<volScalarField> trd
return min
(
new volScalarField
(
min
(
visc
/(
max
(
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)*sqr(kappa_*d)
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
), scalar(10.0)
)
)
visc
/(
max
(
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)*sqr(kappa_*y_)
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
),
scalar(10)
);
return trd;
}
@ -105,43 +102,38 @@ tmp<volScalarField> SpalartAllmarasIDDES::rd
tmp<volScalarField> SpalartAllmarasIDDES::fd(const volScalarField& S) const
{
return 1.0 - tanh(pow3(8.0*rd(nuSgs_+transport().nu(), S)));
return 1 - tanh(pow3(8*rd(nuEff(), S)));
}
void SpalartAllmarasIDDES::dTildaUpdate(const volScalarField& S)
tmp<volScalarField> SpalartAllmarasIDDES::dTilda(const volScalarField& S) const
{
volScalarField alpha = this->alpha();
volScalarField expTerm = exp(sqr(alpha));
volScalarField fHill =
2.0*(pos(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
2*(pos(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
volScalarField fStep = min(2*pow(expTerm, -9.0), scalar(1));
volScalarField fHyb = max(1 - fd(S), fStep);
volScalarField fAmp = 1 - max(ft(S), fl(S));
volScalarField fRestore = max(fHill - 1, scalar(0))*fAmp;
volScalarField fStep = min(2.0*pow(expTerm, -9.0), scalar(1));
volScalarField fHyb = max(1.0 - fd(S), fStep);
volScalarField fAmp = 1.0 - max(ft(S), fl(S));
volScalarField fRestore = max(fHill - 1.0, scalar(0))*fAmp;
// volScalarField ft2 = IGNORING ft2 terms
// IGNORING ft2 terms
volScalarField Psi = sqrt
(
min
(
scalar(100),
(1.0 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*fv2())/max(SMALL, fv1())
(1 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*fv2())/max(SMALL, fv1())
)
);
dTilda_ = max
return max
(
dimensionedScalar("zero", dimLength, 0.0),
fHyb*(1.0 + fRestore*Psi)*wallDist(mesh_).y()
+ (1.0 - fHyb)*CDES_*Psi*delta()
dimensionedScalar("SMALL", dimLength, SMALL),
fHyb*(1 + fRestore*Psi)*y_
+ (1 - fHyb)*CDES_*Psi*delta()
);
}

10
src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.H
View File

@ -69,12 +69,16 @@ class SpalartAllmarasIDDES
tmp<volScalarField> alpha() const;
tmp<volScalarField> ft(const volScalarField& S) const;
tmp<volScalarField> fl(const volScalarField& S) const;
tmp<volScalarField> rd
(
const volScalarField& visc,
const volScalarField& S
) const;
//- Delay function
tmp<volScalarField> fd(const volScalarField& S) const;
// Disallow default bitwise copy construct and assignment
SpalartAllmarasIDDES(const SpalartAllmarasIDDES&);
SpalartAllmarasIDDES& operator=(const SpalartAllmarasIDDES&);
@ -84,11 +88,9 @@ protected:
// Protected member functions
//- Delay function
tmp<volScalarField> fd(const volScalarField& S) const;
//- Length scale
virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
//- Length scale calculation
virtual void dTildaUpdate(const volScalarField& S);
public: