Merge branch 'nextRelease' of ssh://git.code.sf.net/p/foam-extend/foam-extend-3.2 into nextRelease

This commit is contained in:
Hrvoje Jasak 2016-08-01 14:45:02 +01:00
commit 990f036ea4
34 changed files with 846 additions and 205 deletions

7
.gitignore vendored
View file

@ -13,6 +13,9 @@
# file-browser settings - anywhere
.directory
# KDevelop include-path files - anywhere
.kdev_include_paths
# CVS recovered versions - anywhere
.#*
@ -104,8 +107,8 @@ etc/cshrc
etc/bashrc
etc/settings.csh
etc/settings.sh
etc/pref.csh
etc/pref.sh
etc/prefs.csh
etc/prefs.sh
# make sure that this settings file is not used
etc/bashrc.preset

View file

@ -37,9 +37,11 @@ using namespace Foam;
int main(int argc, char *argv[])
{
# include "addRegionOption.H"
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createNamedMesh.H"
# include "createFaMesh.H"
Info<< "Time = " << runTime.timeName() << nl << endl;

View file

@ -50,6 +50,8 @@ using namespace Foam;
int main(int argc, char *argv[])
{
# include "addRegionOption.H"
argList::validArgs.append("STL mesh file");
argList args(argc, argv);

View file

@ -23,11 +23,18 @@
faceRegion[faceI] = patch[faceI].region();
}
word polyRegionName;
if (!args.optionReadIfPresent("region", polyRegionName))
{
polyRegionName = polyMesh::defaultRegion;
}
polyMesh mesh
(
IOobject
(
polyMesh::defaultRegion,
polyRegionName,
runTime.constant(),
runTime
),

View file

@ -63,9 +63,11 @@ public:
int main(int argc, char *argv[])
{
# include "addRegionOption.H"
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createNamedMesh.H"
// Reading faMeshDefinition dictionary
IOdictionary faMeshDefinition

View file

@ -325,7 +325,7 @@ void printAllSets(const polyMesh& mesh, Ostream& os)
// Physically remove a set
void removeSet
bool removeSet
(
const polyMesh& mesh,
const word& setType,
@ -346,7 +346,11 @@ void removeSet
fileName object = objects[setName]->objectPath();
Info<< "Removing file " << object << endl;
rm(object);
return true;
}
return false;
}
@ -392,11 +396,7 @@ bool doCommand
IOobject::readOption r;
if (action == topoSetSource::REMOVE)
{
removeSet(mesh, setType, setName);
}
else if
if
(
(action == topoSetSource::NEW)
|| (action == topoSetSource::CLEAR)
@ -414,7 +414,11 @@ bool doCommand
currentSet.resize(max(currentSet.size(), typSize));
}
if (!currentSetPtr.valid())
if (action == topoSetSource::REMOVE)
{
ok = removeSet(mesh, setType, setName);
}
else if (!currentSetPtr.valid())
{
Info<< " Cannot construct/load set "
<< topoSet::localPath(mesh, setName) << endl;

View file

@ -55,7 +55,7 @@ void faMeshDecomposition::distributeFaces()
(
IOobject
(
fvMesh::defaultRegion,
GeoMesh<polyMesh>::mesh_.name(),
processorDb.timeName(),
processorDb
)
@ -252,7 +252,7 @@ void faMeshDecomposition::decomposeMesh(const bool filterEmptyPatches)
(
IOobject
(
fvMesh::defaultRegion,
GeoMesh<polyMesh>::mesh_.name(),
processorDb.timeName(),
processorDb
)
@ -1073,7 +1073,7 @@ void faMeshDecomposition::decomposeMesh(const bool filterEmptyPatches)
(
IOobject
(
fvMesh::defaultRegion,
GeoMesh<polyMesh>::mesh_.name(),
processorDb.timeName(),
processorDb
)
@ -1191,7 +1191,7 @@ bool faMeshDecomposition::writeDecomposition()
(
IOobject
(
fvMesh::defaultRegion,
GeoMesh<polyMesh>::mesh_.name(),
processorDb.timeName(),
processorDb
)

View file

@ -37,13 +37,29 @@ Foam::autoPtr<Foam::dynamicFvMesh> Foam::dynamicFvMesh::New(const IOobject& io)
forAll(libNames,i) {
const word libName("lib"+libNames[i]+".so");
bool ok=dlLibraryTable::open(libName);
if(!ok) {
WarningIn("dynamicFvMesh::New(const IOobject& io)")
<< "Loading of dynamic mesh library " << libName
bool libLoaded = false;
dlLibraryTable& ll = dlLibraryTable::loadedLibraries;
forAllConstIter(dlLibraryTable, ll, llI)
{
if (ll(llI.key()) == libName)
{
libLoaded = true;
break;
}
}
if (!libLoaded)
{
bool ok=dlLibraryTable::open(libName);
if(!ok) {
WarningIn("dynamicFvMesh::New(const IOobject& io)")
<< "Loading of dynamic mesh library " << libName
<< " unsuccesful. Some dynamic mesh methods may not be "
<< " available"
<< endl;
}
}
}

View file

@ -61,7 +61,7 @@ Foam::UniformDimensionedField<Type>::UniformDimensionedField
{
dictionary dict(readStream(typeName));
this->dimensions().reset(dict.lookup("dimensions"));
this->value() = dict.lookup("value");
this->value() = pTraits<Type>(dict.lookup("value"));
}

View file

@ -425,7 +425,7 @@ public:
//- Return current motion time index
label curMotionTimeIndex() const
{
return moving_;
return curMotionTimeIndex_;
}
//- Set the mesh to be moving

View file

@ -73,6 +73,7 @@ $(faceSources)/boundaryToFace/boundaryToFace.C
$(faceSources)/zoneToFace/zoneToFace.C
$(faceSources)/setToFace/setToFace.C
$(faceSources)/boxToFace/boxToFace.C
$(faceSources)/rotatedBoxToFace/rotatedBoxToFace.C
pointSources = sets/pointSources
$(pointSources)/labelToPoint/labelToPoint.C

View file

@ -46,18 +46,21 @@ addToRunTimeSelectionTable(topoSetSource, pointToCell, istream);
Foam::topoSetSource::addToUsageTable Foam::pointToCell::usage_
(
pointToCell::typeName,
"\n Usage: pointToCell <pointSet> any\n\n"
" Select all cells with any point in the pointSet\n\n"
"\n Usage: pointToCell <pointSet> any|all\n\n"
" Select cells with\n"
" -any point in the pointSet\n"
" -all points in the pointSet\n\n"
);
template<>
const char* Foam::NamedEnum<Foam::pointToCell::pointAction, 1>::names[] =
const char* Foam::NamedEnum<Foam::pointToCell::pointAction, 2>::names[] =
{
"any"
"any",
"all"
};
const Foam::NamedEnum<Foam::pointToCell::pointAction, 1>
const Foam::NamedEnum<Foam::pointToCell::pointAction, 2>
Foam::pointToCell::pointActionNames_;
@ -69,9 +72,9 @@ void Foam::pointToCell::combine(topoSet& set, const bool add) const
pointSet loadedSet(mesh_, setName_);
// Handle any selection
if (option_ == ANY)
{
// Add cells with any point in loadedSet
for
(
pointSet::const_iterator iter = loadedSet.begin();
@ -89,6 +92,53 @@ void Foam::pointToCell::combine(topoSet& set, const bool add) const
}
}
}
else if (option_ == ALL)
{
// Add all cells whose points are all in set.
Map<label> numPoints(loadedSet.size());
forAllConstIter(pointSet, loadedSet, iter)
{
label pointI = iter.key();
const labelList& pCells = mesh_.pointCells()[pointI];
forAll(pCells, pCellI)
{
label cellI = pCells[pCellI];
Map<label>::iterator fndCell = numPoints.find(cellI);
if (fndCell == numPoints.end())
{
numPoints.insert(cellI, 1);
}
else
{
fndCell()++;
}
}
}
// Include cells that are referenced as many times as there are points
// in cell -> all points of cell
for
(
Map<label>::const_iterator iter = numPoints.begin();
iter != numPoints.end();
++iter
)
{
label cellI = iter.key();
if (iter() == mesh_.cellPoints()[cellI].size())
{
addOrDelete(set, cellI, add);
}
}
}
}

View file

@ -55,8 +55,8 @@ public:
//- Enumeration defining the valid options
enum pointAction
{
ANY // Cells using any point in set
//ALL // Possible extension: cells whose all points are in set
ANY, // Cells using any point in set
ALL // Cells whose all points are in set
};
private:
@ -64,7 +64,7 @@ private:
//- Add usage string
static addToUsageTable usage_;
static const NamedEnum<pointAction, 1> pointActionNames_;
static const NamedEnum<pointAction, 2> pointActionNames_;
//- Name of set to use
word setName_;

View file

@ -0,0 +1,193 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "rotatedBoxToFace.H"
#include "polyMesh.H"
#include "cellModeller.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(rotatedBoxToFace, 0);
addToRunTimeSelectionTable(topoSetSource, rotatedBoxToFace, word);
addToRunTimeSelectionTable(topoSetSource, rotatedBoxToFace, istream);
}
Foam::topoSetSource::addToUsageTable Foam::rotatedBoxToFace::usage_
(
rotatedBoxToFace::typeName,
"\n Usage: rotatedBoxToFace (originx originy originz)"
" (ix iy iz) (jx jy jz) (kx ky kz)\n\n"
" Select all faces with faceCentre within parallelopiped\n\n"
);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::rotatedBoxToFace::combine(topoSet& set, const bool add) const
{
// Define a cell for the box
pointField boxPoints(8);
boxPoints[0] = origin_;
boxPoints[1] = origin_ + i_;
boxPoints[2] = origin_ + i_ + j_;
boxPoints[3] = origin_ + j_;
boxPoints[4] = origin_ + k_;
boxPoints[5] = origin_ + k_ + i_;
boxPoints[6] = origin_ + k_ + i_ + j_;
boxPoints[7] = origin_ + k_ + j_;
labelList boxVerts(8);
forAll(boxVerts, i)
{
boxVerts[i] = i;
}
const cellModel& hex = *(cellModeller::lookup("hex"));
// Get outwards pointing faces.
faceList boxFaces(cellShape(hex, boxVerts).faces());
// Precalculate normals
vectorField boxFaceNormals(boxFaces.size());
forAll(boxFaces, i)
{
boxFaceNormals[i] = boxFaces[i].normal(boxPoints);
Pout<< "Face:" << i << " position:" << boxFaces[i].centre(boxPoints)
<< " normal:" << boxFaceNormals[i] << endl;
}
// Check whether face centre is inside all faces of box.
const pointField& ctrs = mesh_.faceCentres();
forAll(ctrs, faceI)
{
bool inside = true;
forAll(boxFaces, i)
{
const face& f = boxFaces[i];
if (((ctrs[faceI] - boxPoints[f[0]]) & boxFaceNormals[i]) > 0)
{
inside = false;
break;
}
}
if (inside)
{
addOrDelete(set, faceI, add);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::rotatedBoxToFace::rotatedBoxToFace
(
const polyMesh& mesh,
const vector& origin,
const vector& i,
const vector& j,
const vector& k
)
:
topoSetSource(mesh),
origin_(origin),
i_(i),
j_(j),
k_(k)
{}
// Construct from dictionary
Foam::rotatedBoxToFace::rotatedBoxToFace
(
const polyMesh& mesh,
const dictionary& dict
)
:
topoSetSource(mesh),
origin_(dict.lookup("origin")),
i_(dict.lookup("i")),
j_(dict.lookup("j")),
k_(dict.lookup("k"))
{}
// Construct from Istream
Foam::rotatedBoxToFace::rotatedBoxToFace(const polyMesh& mesh, Istream& is)
:
topoSetSource(mesh),
origin_(is),
i_(is),
j_(is),
k_(is)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::rotatedBoxToFace::~rotatedBoxToFace()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::rotatedBoxToFace::applyToSet
(
const topoSetSource::setAction action,
topoSet& set
) const
{
if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
{
Info<< " Adding faces with center within rotated box " << endl;
combine(set, true);
}
else if (action == topoSetSource::DELETE)
{
Info<< " Removing faces with center within rotated box " << endl;
combine(set, false);
}
}
// ************************************************************************* //

View file

@ -0,0 +1,133 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::rotatedBoxToFace
Description
A topoSetSource to select faces based on face centres inside
rotated/skewed box (parallelopiped?).
Box defined as origin and i,j,k vectors.
E.g. box rotated 45 degrees around z-axis with sizes sqrt(0.2^2+0.2^2)
(and extra large, 200 in z direction):
@verbatim
origin ( 0.4 0.4 -100);
i ( 0.2 0.2 0);
j (-0.2 0.2 0);
k ( 0.0 0.0 100);
@endverbatim
SourceFiles
rotatedBoxToFace.C
\*---------------------------------------------------------------------------*/
#ifndef rotatedBoxToFace_H
#define rotatedBoxToFace_H
#include "topoSetSource.H"
#include "treeBoundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class rotatedBoxToFace Declaration
\*---------------------------------------------------------------------------*/
class rotatedBoxToFace
:
public topoSetSource
{
// Private data
//- Add usage string
static addToUsageTable usage_;
//- skewed box
const vector origin_;
const vector i_;
const vector j_;
const vector k_;
// Private Member Functions
void combine(topoSet& set, const bool add) const;
public:
//- Runtime type information
TypeName("rotatedBoxToFace");
// Constructors
//- Construct from components
rotatedBoxToFace
(
const polyMesh& mesh,
const vector& origin,
const vector& i,
const vector& j,
const vector& k
);
//- Construct from dictionary
rotatedBoxToFace(const polyMesh& mesh, const dictionary& dict);
//- Construct from Istream
rotatedBoxToFace(const polyMesh& mesh, Istream&);
// Destructor
virtual ~rotatedBoxToFace();
// Member Functions
virtual void applyToSet
(
const topoSetSource::setAction action,
topoSet&
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -17,6 +17,7 @@ cyclicGgiTetPolyPatch = $(constraintTetPolyPatches)/cyclicGgi
mixingPlaneTetPolyPatch = $(constraintTetPolyPatches)/mixingPlane
globalTetPolyPatch = $(constraintTetPolyPatches)/global
wallTetPolyPatch = $(derivedTetPolyPatches)/wall
directMappedTetPolyPatch = $(derivedTetPolyPatches)/directMapped
MapTetFemFields = $(tetPolyMesh)/MapTetFemFields
@ -39,6 +40,8 @@ $(mixingPlaneTetPolyPatch)/mixingPlaneTetPolyPatch.C
$(globalTetPolyPatch)/globalTetPolyPatch.C
$(globalTetPolyPatch)/calcGlobalTetPolyPatchAddr.C
$(wallTetPolyPatch)/wallTetPolyPatch.C
$(directMappedTetPolyPatch)/directMappedTetPolyPatch.C
$(directMappedTetPolyPatch)/directMappedWallTetPolyPatch.C
$(tetPolyBoundaryMesh)/tetPolyBoundaryMesh.C
$(tetPolyMesh)/tetPolyMeshLduAddressing.C
$(tetPolyMesh)/tetPolyMesh.C

View file

@ -0,0 +1,52 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "directMappedTetPolyPatch.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
defineTypeNameAndDebug(directMappedTetPolyPatch, 0);
// Add the patch constructor functions to the hash tables
addToRunTimeSelectionTable
(
faceTetPolyPatch,
directMappedTetPolyPatch,
polyPatch
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,86 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::directMappedTetPolyPatch
Description
Direct mapped patch.
SourceFiles
directMappedTetPolyPatch.C
\*---------------------------------------------------------------------------*/
#ifndef directMappedTetPolyPatch_H
#define directMappedTetPolyPatch_H
#include "faceTetPolyPatch.H"
#include "directMappedPolyPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class directMappedTetPolyPatch Declaration
\*---------------------------------------------------------------------------*/
class directMappedTetPolyPatch
:
public faceTetPolyPatch
{
public:
//- Runtime type information
TypeName(directMappedPolyPatch::typeName_());
// Constructors
//- Construct from polyPatch
directMappedTetPolyPatch
(
const polyPatch& patch,
const tetPolyBoundaryMesh& bm
)
:
faceTetPolyPatch(patch, bm)
{}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,52 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "directMappedWallTetPolyPatch.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
defineTypeNameAndDebug(directMappedWallTetPolyPatch, 0);
// Add the patch constructor functions to the hash tables
addToRunTimeSelectionTable
(
faceTetPolyPatch,
directMappedWallTetPolyPatch,
polyPatch
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,84 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::directMappedWallTetPolyPatch
Description
Direct mapped wall patch.
SourceFiles
directMappedWallTetPolyPatch.C
\*---------------------------------------------------------------------------*/
#ifndef directMappedWallTetPolyPatch_H
#define directMappedWallTetPolyPatch_H
#include "wallTetPolyPatch.H"
#include "directMappedWallPolyPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class directMappedWallTetPolyPatch Declaration
\*---------------------------------------------------------------------------*/
class directMappedWallTetPolyPatch
:
public wallTetPolyPatch
{
public:
//- Runtime type information
TypeName(directMappedWallPolyPatch::typeName_());
// Constructors
//- Construct from polyPatch
directMappedWallTetPolyPatch
(
const polyPatch& patch,
const tetPolyBoundaryMesh& bm
)
:
wallTetPolyPatch(patch, bm)
{}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -56,12 +56,12 @@ Foam::aungierRedlichKwong::aungierRedlichKwong(Istream& is)
//CL: Only uses the default values
rhoMin_(1e-3),
rhoMax_(1500),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
TSave(0.0),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
{
is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
}
@ -91,12 +91,12 @@ Foam::aungierRedlichKwong::aungierRedlichKwong(const dictionary& dict)
//CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
{
is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
}

View file

@ -87,9 +87,6 @@ class aungierRedlichKwong
scalar rhoMin_;
scalar rhoMax_;
//- Density @STD, initialise after a0, b, c, n!
scalar rhostd_;
//CL: Variables to save the values of a, dadT and d2adT2 of the mixture
//CL: Variables must corrected for changing temperatures
mutable scalar aSave;
@ -99,6 +96,9 @@ class aungierRedlichKwong
//CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct
mutable scalar TSave;
//- Density @STD, initialise after a0, b, c, n!
scalar rhostd_;
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const;

View file

@ -68,11 +68,11 @@ inline aungierRedlichKwong::aungierRedlichKwong(const word& name, const aungierR
c2_(pg.c2_),
rhoMin_(pg.rhoMin_),
rhoMax_(pg.rhoMax_),
rhostd_(pg.rhostd_),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
TSave(0.0),
rhostd_(pg.rhostd_)
{}
@ -237,13 +237,16 @@ inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
scalar cpVm = c_ + Vm;
return
(
a(T)*(b3_ - 2*b2_*c_ + b_*(c_ + Vm)*(c_ - 3*Vm) + 2*Vm*pow(c_ + Vm, 2))
- this->RR()*T*Vm2*(b2_ + 2*b_*Vm + Vm2)
a(T)*(b3_ - 2*b2_*c_ + b_*cpVm*(c_ - 3*Vm) + 2*Vm*cpVm*cpVm)
- this->RR()*T*Vm2*bpVm2
)
/(Vm2*pow(b_ - c_ - Vm, 2)*pow(b_ + Vm, 2));
/(Vm2*pow(b_ - c_ - Vm, 2)*bpVm2);
}
@ -284,8 +287,7 @@ inline scalar aungierRedlichKwong::integral_p_dv
{
scalar Vm = this->W()/rho;
return this->RR()*T*log(-b_ + c_ + Vm) + a(T)/b_*(log(b_ + Vm) - log(Vm));
//return this->RR()*T*log(-b_ + c_ + Vm) + a(T)*log(b_ + Vm)/b_ - a(T)*log(Vm)/b_;
return this->RR()*T*log(Vm - b_ + c_) + a(T)/b_*log((b_ + Vm)/Vm);
}
@ -299,8 +301,7 @@ inline scalar aungierRedlichKwong::integral_dpdT_dv
{
scalar Vm = this->W()/rho;
return this->RR()*log(-b_ + c_ + Vm) + dadT(T)/b_*(log(b_ + Vm) - log(Vm));
//return this->RR()*log(-b_ + c_ + Vm) + dadT(T)*log(b_ + Vm)/b_ - dadT(T)*log(Vm)/b_;
return this->RR()*log(Vm - b_ + c_) + dadT(T)/b_*log((b_ + Vm)/Vm);
}
@ -309,7 +310,7 @@ inline scalar aungierRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
{
scalar Vm = this->W()/rho;
return -d2adT2(T)/(Vm*(Vm + b_));
return -d2adT2(T)/(Vm*(b_ + Vm));
}
@ -400,7 +401,7 @@ inline scalar aungierRedlichKwong::integral_d2pdT2_dv
{
scalar Vm = this->W()/rho;
return d2adT2(T)*log(b_ + Vm)/b_ - d2adT2(T)*log(Vm)/b_;
return d2adT2(T)/b_*log((b_ + Vm)/Vm);
}

View file

@ -47,19 +47,15 @@ Foam::pengRobinson::pengRobinson(Istream& is)
b_(0.077796*this->RR()*Tcrit_/pcrit_),
n_(0.37464 + 1.54226*azentricFactor_ - 0.26992*pow(azentricFactor_, 2)),
b2_(b_*b_),
b3_(b2_*b_),
b4_(b3_*b_),
b5_(b4_*b_),
b6_(b5_*b_),
//CL: Only uses the default values
rhoMin_(1e-3),
rhoMax_(1500),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
TSave(0.0),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
{
is.check("pengRobinson::pengRobinson(Istream& is)");
}

View file

@ -79,18 +79,11 @@ class pengRobinson
//CL: pow of constants b_ used in the code e.g. b2_=b*b;
scalar b2_;
scalar b3_;
scalar b4_;
scalar b5_;
scalar b6_;
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
scalar rhoMin_;
scalar rhoMax_;
//- Density @STD, initialise after a0, b!
scalar rhostd_;
//CL: Variables to save the values of a, dadT and d2adT2 of the mixture
//CL: Variables must corrected for changing temperatures
mutable scalar aSave;
@ -100,6 +93,9 @@ class pengRobinson
//CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct
mutable scalar TSave;
//- Density @STD, initialise after a0, b!
scalar rhostd_;
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const;

View file

@ -60,17 +60,13 @@ inline pengRobinson::pengRobinson(const word& name, const pengRobinson& pr)
b_(pr.b_),
n_(pr.n_),
b2_(pr.b2_),
b3_(pr.b3_),
b4_(pr.b4_),
b5_(pr.b5_),
b6_(pr.b6_),
rhoMin_(pr.rhoMin_),
rhoMax_(pr.rhoMax_),
rhostd_(pr.rhostd_),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
TSave(0.0),
rhostd_(pr.rhostd_)
{}
@ -123,7 +119,7 @@ inline void pengRobinson::updateModelCoefficients(const scalar T) const
d2aSave = a0_*n_*(n_ + 1)*TTc/(2*T*T);
//CL: saving the temperature at which the coefficients are valid
TSave=T;
TSave = T;
}
@ -131,7 +127,7 @@ inline void pengRobinson::updateModelCoefficients(const scalar T) const
inline scalar pengRobinson::a(const scalar T) const
{
//CL: check if a has already been calculated for this temperature
if(TSave==T)
if(TSave == T)
{
return aSave;
}
@ -148,7 +144,7 @@ inline scalar pengRobinson::a(const scalar T) const
inline scalar pengRobinson::dadT(const scalar T) const
{
// check if a has already been calculated for this temperature
if(TSave==T)
if(TSave == T)
{
return daSave;
}
@ -165,7 +161,7 @@ inline scalar pengRobinson::dadT(const scalar T) const
inline scalar pengRobinson::d2adT2(const scalar T) const
{
// check if a has already been calculated for this temperature
if(TSave==T)
if(TSave == T)
{
return d2aSave;
}
@ -210,22 +206,15 @@ inline scalar pengRobinson::p(const scalar rho, const scalar T) const
inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
return
(
2*a(T)*
(
b3_ - b2_*Vm - b_*Vm2 + Vm3
)
- this->RR()*T*
(
b4_ - 4*b3_*Vm + 2*b2_*Vm2 + 4*b_*Vm3 + Vm4
)
)
/pow(b_ - Vm, 8);
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bpVm = b_ + Vm;
scalar t1 = bmVm2 - 2*Vm*Vm;
scalar t12 = t1*t1;
return (2*a(T)*bmVm2*bpVm - this->RR()*T*t12)/(bmVm2*t12);
}
@ -285,7 +274,7 @@ inline scalar pengRobinson::integral_dpdT_dv
) const
{
scalar Vm = this->W()/rho;
scalar root2=pow(2, 0.5);
scalar root2 = pow(2, 0.5);
return
- root2*dadT(T)*log(b_*(1 - root2) + Vm)/(4*b_)
@ -308,23 +297,16 @@ inline scalar pengRobinson::d2pdv2(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
scalar Vm5 = Vm4*Vm;
scalar Vm6 = Vm5*Vm;
return 2*
(
a(T)*
(
5*b5_ - 9*b4_*Vm + 4*b2_*Vm3 + 3*b_*Vm4 - 3*Vm5
)
- this->RR()*T*
(
b6_ - 6*b5_*Vm + 9*b4_*Vm2 + 4*b3_*Vm3 - 9*b2_*Vm4 - 6*b_*Vm5 - Vm6
)
)
/pow(b_ - Vm, 9);
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bmVm3 = bmVm2*bmVm;
scalar t1 = bmVm2 - 2*Vm2;
scalar t12 = t1*t1;
scalar t13 = t12*t1;
return 2*(a(T)*bmVm3*(5*b_ + 6*b_*Vm + 3*Vm2) - this->RR()*T*t13)/(bmVm*t1);
}
@ -360,16 +342,16 @@ inline scalar pengRobinson::d2pdvdT
) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
return
(
2*dadT(T)*(b3_ - b2_*Vm - b_*Vm2 + Vm3)
- this->RR()*(b4_ - 4*b3_*Vm + 2*b2_*Vm2 + 4*b_*Vm3 + Vm4)
)
/pow(b_ - Vm, 6);
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bpVm = b_ + Vm;
scalar t1 = bmVm2 - 2*Vm*Vm;
scalar t12 = t1*t1;
scalar t14 = t12*t12;
return (2*dadT(T)*bmVm2*bpVm - this->RR()*t14)/(bmVm2*t12);
}

View file

@ -45,8 +45,6 @@ Foam::redlichKwong::redlichKwong(Istream& is)
a_(0.42748*pow(this->RR(), 2)*pow(Tcrit_, 2.5)/pcrit_),
b_(0.08664*this->RR()*Tcrit_/pcrit_),
b2_(pow(b_,2)),
b3_(pow(b_,3)),
b5_(pow(b_,5)),
//CL: Only uses the default values
rhoMin_(1e-3),
rhoMax_(1500),

View file

@ -69,8 +69,6 @@ class redlichKwong
//CL: pow of constants b_ used in the code e.g. b2_=b*b;
scalar b2_;
scalar b3_;
scalar b5_;
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
scalar rhoMin_;

View file

@ -57,8 +57,6 @@ inline redlichKwong::redlichKwong(const word& name, const redlichKwong& rk)
a_(rk.a_),
b_(rk.b_),
b2_(rk.b2_),
b3_(rk.b3_),
b5_(rk.b5_),
rhoMin_(rk.rhoMin_),
rhoMax_(rk.rhoMax_),
rhostd_(rk.rhostd_)
@ -132,15 +130,13 @@ inline scalar redlichKwong::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar T05 = pow(T, 0.5);
scalar bpVm2 = pow(b_ + Vm, 2);
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
return
(
a_*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*T*T05*Vm2*bpVm2
)
/(T05*Vm2*bpVm2*pow(b_ - Vm, 2));
return (a_*bmVm2*bpVm - this->RR()*T*T05*Vm2*bpVm2)/(T05*Vm2*bpVm2*bmVm2);
}
@ -181,8 +177,7 @@ inline scalar redlichKwong::integral_p_dv
{
scalar Vm = this->W()/rho;
return this->RR()*T*log(Vm - b_) + a_/(b_*sqrt(T))*(log(b_ + Vm) - log(Vm));
//return this->RR()*T*log(Vm - b_) + (a_*log(b_ + Vm))/(b_*sqrt(T)) - (a_*log(Vm))/(b_*sqrt(T));
return this->RR()*T*log(Vm - b_) + a_/(b_*sqrt(T))*log((b_ + Vm)/Vm);
}
@ -196,8 +191,7 @@ inline scalar redlichKwong::integral_dpdT_dv
{
scalar Vm = this->W()/rho;
return this->RR()*log(Vm - b_) - a_/(2*b_*pow(T, 1.5))*(log(b_ + Vm) - log(Vm));
//return this->RR()*log(Vm - b_) - (a_*log(b_ + Vm))/(2*b_*T15_) + (a_*log(Vm))/(2*b_*T15_);
return this->RR()*log(Vm - b_) - a_/(2*b_*pow(T, 1.5))*log((b_ + Vm)/Vm);
}
@ -205,6 +199,7 @@ inline scalar redlichKwong::integral_dpdT_dv
inline scalar redlichKwong::d2pdT2(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
return -0.75*a_/(pow(T, 2.5)*Vm*(b_ + Vm));
}
@ -215,19 +210,16 @@ inline scalar redlichKwong::d2pdv2(const scalar rho, const scalar T) const
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
scalar Vm5 = Vm4*Vm;
scalar T05 = pow(T, 0.5);
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bmVm3 = bmVm2*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
scalar bpVm3 = bpVm2*bpVm;
return
(
2*
(
a_*(b5_ - 3*b3_*Vm2 - b2_*Vm3 + 6*b_*Vm4 - 3*Vm5)
+ this->RR()*T*T05*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3)
)
/(T05*Vm3*pow((b_ + Vm), 3)*pow(Vm - b_, 3))
);
return 2*(a_*bmVm3*(b2_ + 3*b_*Vm + 3*Vm2) + this->RR()*T*T05*Vm3*bpVm3)
/(T05*Vm3*bpVm3*bmVm3);
}
@ -256,18 +248,14 @@ inline scalar redlichKwong::d2pdvdT(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar T15_ = pow(T, 1.5);
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
return
-(
0.5*
(
a_*(b3_ - 3*b_*Vm2 + 2*Vm3)
+ 2*this->RR()*T15_*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
)
/(T15_*Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
return -(0.5*a_*bmVm2*(b_ + 2*Vm) + this->RR()*T15_*Vm2*bpVm2)
/(T15_*Vm2*bpVm2*bmVm2);
}
@ -279,9 +267,10 @@ inline scalar redlichKwong::integral_d2pdT2_dv
const scalar T
) const
{
scalar T25_=pow(T,2.5);
scalar Vm = this->W()/rho;
return 0.75*a_*log(b_ + Vm)/(T25_*b_) - 0.75*a_*log(Vm)/(T25_*b_);
scalar T25 = pow(T, 2.5);
return 0.75*a_/(T25*b_)*log((b_ + Vm)/Vm);
}

View file

@ -47,17 +47,15 @@ Foam::soaveRedlichKwong::soaveRedlichKwong(Istream& is)
b_(0.08664*this->RR()*Tcrit_/pcrit_),
n_(0.48508 + 1.55171*azentricFactor_ - 0.15613*pow(azentricFactor_, 2)),
b2_(b_*b_),
b3_(b2_*b_),
b5_(b2_*b3_),
//CL: Only uses the default values
rhoMin_(1e-3),
rhoMax_(1500),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
TSave(0.0),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
{
is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)");
}
@ -83,12 +81,12 @@ soaveRedlichKwong::soaveRedlichKwong(const dictionary& dict)
//CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
TSave(0.0),
// Starting GUESS for the density by ideal gas law
rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
{
is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)");
}

View file

@ -76,16 +76,11 @@ class soaveRedlichKwong
//CL: pow of constants b_ used in the code e.g. b2_=b*b;
scalar b2_;
scalar b3_;
scalar b5_;
//CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
scalar rhoMin_;
scalar rhoMax_;
//- Density @STD, initialise after a0, b!
scalar rhostd_;
//CL: Variables to save the values of a, dadT and d2adT2 of the mixture
//CL: Variables must corrected for changing temperatures
mutable scalar aSave;
@ -95,6 +90,9 @@ class soaveRedlichKwong
//CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct
mutable scalar TSave;
//- Density @STD, initialise after a0, b!
scalar rhostd_;
//Protected functions
//CL: function updates the coefficients (aSave, daSave, d2aSave)
inline void updateModelCoefficients(const scalar T) const;

View file

@ -60,15 +60,13 @@ inline soaveRedlichKwong::soaveRedlichKwong(const word& name, const soaveRedlich
b_(srk.b_),
n_(srk.n_),
b2_(srk.b2_),
b3_(srk.b3_),
b5_(srk.b5_),
rhoMin_(srk.rhoMin_),
rhoMax_(srk.rhoMax_),
rhostd_(srk.rhostd_),
aSave(0.0),
daSave(0.0),
d2aSave(0.0),
TSave(0.0)
TSave(0.0),
rhostd_(srk.rhostd_)
{}
@ -129,7 +127,7 @@ inline void soaveRedlichKwong::updateModelCoefficients(const scalar T) const
inline scalar soaveRedlichKwong::a(const scalar T) const
{
//CL: check if a has already been calculated for this temperature
if(TSave==T)
if(TSave == T)
{
return aSave;
}
@ -146,7 +144,7 @@ inline scalar soaveRedlichKwong::a(const scalar T) const
inline scalar soaveRedlichKwong::dadT(const scalar T) const
{
// check if a has already been calculated for this temperature
if(TSave==T)
if(TSave == T)
{
return daSave;
}
@ -163,7 +161,7 @@ inline scalar soaveRedlichKwong::dadT(const scalar T) const
inline scalar soaveRedlichKwong::d2adT2(const scalar T) const
{
// check if a has already been calculated for this temperature
if(TSave==T)
if(TSave == T)
{
return d2aSave;
}
@ -208,13 +206,12 @@ inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
return
(
a(T)*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*T*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
/(Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
return(a(T)*bmVm2*(b_ + 2*Vm) - this->RR()*T*Vm2*bpVm2)/(Vm2*bpVm2*bmVm2);
}
@ -255,8 +252,7 @@ inline scalar soaveRedlichKwong::integral_p_dv
{
scalar Vm = this->W()/rho;
return this->RR()*T*log(Vm - b_) + a(T)/b_*(log(b_ + Vm) - log(Vm));
//return this->RR()*T*log(Vm - b_) + a(T)*log(b_ + Vm)/b_ - a(T)*log(Vm)/b_;
return this->RR()*T*log(Vm - b_) + a(T)/b_*log((b_ + Vm)/Vm);
}
@ -271,8 +267,7 @@ inline scalar soaveRedlichKwong::integral_dpdT_dv
{
scalar Vm = this->W()/rho;
return this->RR()*log(Vm - b_) + dadT(T)/b_*(log(b_ + Vm) - log(Vm));
//return this->RR()*log(Vm - b_) + dadT(T)*log(b_ + Vm)/b_ - dadT(T)*log(Vm)/b_;
return this->RR()*log(Vm - b_) + dadT(T)/b_*log((b_ + Vm)/Vm);
}
@ -291,15 +286,15 @@ inline scalar soaveRedlichKwong::d2pdv2(const scalar rho, const scalar T) const
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar Vm4 = Vm3*Vm;
scalar Vm5 = Vm4*Vm;
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bmVm3 = bmVm2*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
scalar bpVm3 = bpVm2*bpVm;
return 2*
(
a(T)*(b5_ - 3*b3_*Vm2 - b2_*Vm3 + 6*b_*Vm4 - 3*Vm5)
+ this->RR()*T*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3)
)
/(Vm3*pow(b_ + Vm, 3)*pow(Vm - b_, 3));
return 2*(a(T)*bmVm3*(b2_ + 3*b_*Vm + 3*Vm2) + this->RR()*T*Vm3*bpVm3)
/(Vm3*bpVm3*bmVm3);
}
@ -336,13 +331,13 @@ inline scalar soaveRedlichKwong::d2pdvdT
{
scalar Vm = this->W()/rho;
scalar Vm2 = Vm*Vm;
scalar Vm3 = Vm2*Vm;
scalar bmVm = b_ - Vm;
scalar bmVm2 = bmVm*bmVm;
scalar bpVm = b_ + Vm;
scalar bpVm2 = bpVm*bpVm;
return
(
dadT(T)*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*Vm2*(b2_ + 2*b_*Vm + Vm2)
)
/(Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
return (dadT(T)*bmVm2*(b_ + 2*Vm) - this->RR()*Vm2*bmVm2)
/(Vm2*bpVm2*bmVm2);
}
@ -356,7 +351,7 @@ inline scalar soaveRedlichKwong::integral_d2pdT2_dv
{
scalar Vm = this->W()/rho;
return d2adT2(T)*log(b_ + Vm)/b_ - d2adT2(T)*log(Vm)/b_;
return d2adT2(T)/b_*log((b_ + Vm)/Vm);
}

View file

@ -18,14 +18,14 @@ FoamFile
//CL: List of possible real gas models
thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 467.6 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801E-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
//thermoType realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
//mixture CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;

View file

@ -45,7 +45,7 @@ timePrecision 10;
adjustTimeStep yes;
maxCo 0.25;
maxCo 0.1;
maxDeltaT 1e-2;