Merge branch 'master' of github.com:OpenFOAM/OpenFOAM-2.3.x

This commit is contained in:
sergio 2014-07-10 11:01:20 +01:00
commit c15e69ea88
27 changed files with 449 additions and 119 deletions

View file

@ -40,8 +40,6 @@ if (pimple.transonic())
fvOptions(psi, p, rho.name())
);
fvOptions.constrain(pEqn);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
@ -74,8 +72,6 @@ else
fvOptions(psi, p, rho.name())
);
fvOptions.constrain(pEqn);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -24,6 +24,9 @@ License
Application
rhoPimpleFoam
Group
grpCompressibleSolvers grpMovingMeshSolvers
Description
Transient solver for laminar or turbulent flow of compressible fluids
for HVAC and similar applications.
@ -71,17 +74,21 @@ int main(int argc, char *argv[])
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
{
// Store divrhoU from the previous time-step/mesh for the correctPhi
volScalarField divrhoU
(
"divrhoU",
fvc::div(fvc::absolute(phi, rho, U))
);
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// Store momentum to set rhoUf for introduced faces.
volVectorField rhoU("rhoU", rho*U);
// Store divrhoU from the previous time-step/mesh for the correctPhi
volScalarField divrhoU(fvc::div(fvc::absolute(phi, rho, U)));
// Do any mesh changes
mesh.update();
@ -102,12 +109,9 @@ int main(int argc, char *argv[])
#include "meshCourantNo.H"
}
if (pimple.nCorrPIMPLE() <= 1)
{
#include "rhoEqn.H"
Info<< "rhoEqn max/min : " << max(rho).value()
<< " " << min(rho).value() << endl;
}
#include "rhoEqn.H"
Info<< "rhoEqn max/min : " << max(rho).value()
<< " " << min(rho).value() << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -47,6 +47,7 @@ Description
#include "fvMesh.H"
#include "MeshedSurfaces.H"
#include "globalIndex.H"
#include "cellSet.H"
#include "extrudedMesh.H"
#include "extrudeModel.H"
@ -192,6 +193,23 @@ void updateFaceLabels(const mapPolyMesh& map, labelList& faceLabels)
}
void updateCellSet(const mapPolyMesh& map, labelHashSet& cellLabels)
{
const labelList& reverseMap = map.reverseCellMap();
labelHashSet newCellLabels(2*cellLabels.size());
forAll(cellLabels, i)
{
label oldCellI = cellLabels[i];
if (reverseMap[oldCellI] >= 0)
{
newCellLabels.insert(reverseMap[oldCellI]);
}
}
cellLabels.transfer(newCellLabels);
}
int main(int argc, char *argv[])
@ -270,6 +288,9 @@ int main(int argc, char *argv[])
word backPatchName;
labelList backPatchFaces;
// Optional added cells (get written to cellSet)
labelHashSet addedCellsSet;
if (mode == PATCH || mode == MESH)
{
if (flipNormals && mode == MESH)
@ -670,11 +691,33 @@ int main(int argc, char *argv[])
map().reverseFaceMap(),
backPatchFaces
);
// Store added cells
if (mode == MESH)
{
const labelListList addedCells
(
layerExtrude.addedCells
(
meshFromMesh,
layerExtrude.layerFaces()
)
);
forAll(addedCells, faceI)
{
const labelList& aCells = addedCells[faceI];
forAll(aCells, i)
{
addedCellsSet.insert(aCells[i]);
}
}
}
}
else
{
// Read from surface
fileName surfName(dict.lookup("surface"));
surfName.expand();
Info<< "Extruding surfaceMesh read from file " << surfName << nl
<< endl;
@ -810,6 +853,7 @@ int main(int argc, char *argv[])
// Update stored data
updateFaceLabels(map(), frontPatchFaces);
updateFaceLabels(map(), backPatchFaces);
updateCellSet(map(), addedCellsSet);
// Move mesh (if inflation used)
if (map().hasMotionPoints())
@ -913,6 +957,9 @@ int main(int argc, char *argv[])
// Update fields
mesh.updateMesh(map);
// Update local data
updateCellSet(map(), addedCellsSet);
// Move mesh (if inflation used)
if (map().hasMotionPoints())
{
@ -929,6 +976,21 @@ int main(int argc, char *argv[])
<< exit(FatalError);
}
// Need writing cellSet
label nAdded = returnReduce(addedCellsSet.size(), sumOp<label>());
if (nAdded > 0)
{
cellSet addedCells(mesh, "addedCells", addedCellsSet);
Info<< "Writing added cells to cellSet " << addedCells.name()
<< nl << endl;
if (!addedCells.write())
{
FatalErrorIn(args.executable()) << "Failed writing cellSet"
<< addedCells.name()
<< exit(FatalError);
}
}
Info<< "End\n" << endl;
return 0;

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -128,6 +128,21 @@ public:
//- Type of values the UList contains.
typedef T value_type;
//- Type that can be used for storing into
// UList::value_type objects.
typedef T& reference;
//- Type that can be used for storing into
// constant UList::value_type objects
typedef const T& const_reference;
//- The type that can represent the difference between any two
// UList iterator objects.
typedef label difference_type;
//- The type that can represent the size of a UList.
typedef label size_type;
// Ostream operator

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -61,8 +61,8 @@ void Foam::polyPatch::movePoints(PstreamBuffers&, const pointField& p)
void Foam::polyPatch::updateMesh(PstreamBuffers&)
{
primitivePatch::clearGeom();
clearAddressing();
primitivePatch::clearOut();
}
@ -348,6 +348,8 @@ const Foam::labelList& Foam::polyPatch::meshEdges() const
void Foam::polyPatch::clearAddressing()
{
primitivePatch::clearTopology();
primitivePatch::clearPatchMeshAddr();
deleteDemandDrivenData(faceCellsPtr_);
deleteDemandDrivenData(mePtr_);
}

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -92,6 +92,7 @@ clearTopology()
deleteDemandDrivenData(pointEdgesPtr_);
deleteDemandDrivenData(pointFacesPtr_);
deleteDemandDrivenData(edgeLoopsPtr_);
deleteDemandDrivenData(localPointOrderPtr_);
}

View file

@ -150,7 +150,7 @@ Foam::vector Foam::eigenValues(const tensor& t)
if (R2 < Q3)
{
scalar sqrtQ = sqrt(Q);
scalar theta = acos(R/(Q*sqrtQ));
scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ))));
scalar m2SqrtQ = -2*sqrtQ;
scalar aBy3 = a/3;
@ -344,7 +344,7 @@ Foam::vector Foam::eigenValues(const symmTensor& t)
if (R2 < Q3)
{
scalar sqrtQ = sqrt(Q);
scalar theta = acos(R/(Q*sqrtQ));
scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ))));
scalar m2SqrtQ = -2*sqrtQ;
scalar aBy3 = a/3;

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -286,7 +286,7 @@ void Foam::reduce
Pout<< "UPstream::allocateRequest for non-blocking reduce"
<< " : request:" << requestID
<< endl;
}
#else
// Non-blocking not yet implemented in mpi
reduce(Value, bop, tag, communicator);

View file

@ -843,7 +843,9 @@ tmp<surfaceScalarField> CoEulerDdtScheme<Type>::meshPhi
(
"meshPhi",
mesh().time().timeName(),
mesh()
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh(),
dimensionedScalar("0", dimVolume/dimTime, 0.0)

View file

@ -978,7 +978,7 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
fvm.source() =
(
rDtCoef*rho.internalField()*vf.oldTime().internalField()
rDtCoef*rho.oldTime().internalField()*vf.oldTime().internalField()
+ offCentre_(ddt0.internalField())
)*mesh().V0();
}
@ -1079,8 +1079,8 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
fvm.source() =
(
rDtCoef
*alpha.internalField()
*rho.internalField()
*alpha.oldTime().internalField()
*rho.oldTime().internalField()
*vf.oldTime().internalField()
+ offCentre_(ddt0.internalField())
)*mesh().V0();
@ -1426,7 +1426,22 @@ tmp<surfaceScalarField> CrankNicolsonDdtScheme<Type>::meshPhi
coef0_(meshPhi0)*mesh().phi().oldTime() - offCentre_(meshPhi0());
}
return coef_(meshPhi0)*mesh().phi() - offCentre_(meshPhi0());
return tmp<surfaceScalarField>
(
new surfaceScalarField
(
IOobject
(
mesh().phi().name(),
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
coef_(meshPhi0)*mesh().phi() - offCentre_(meshPhi0())
)
);
}

View file

@ -845,7 +845,10 @@ tmp<surfaceScalarField> SLTSDdtScheme<Type>::meshPhi
(
"meshPhi",
mesh().time().timeName(),
mesh()
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh(),
dimensionedScalar("0", dimVolume/dimTime, 0.0)

View file

@ -976,7 +976,22 @@ tmp<surfaceScalarField> backwardDdtScheme<Type>::meshPhi
// Coefficient for t-1/2 (between times n and 0)
scalar coefftn_0 = 1 + coefft0_00;
return coefftn_0*mesh().phi() - coefft0_00*mesh().phi().oldTime();
return tmp<surfaceScalarField>
(
new surfaceScalarField
(
IOobject
(
mesh().phi().name(),
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
coefftn_0*mesh().phi() - coefft0_00*mesh().phi().oldTime()
)
);
}

View file

@ -740,7 +740,10 @@ tmp<surfaceScalarField> localEulerDdtScheme<Type>::meshPhi
(
"meshPhi",
mesh().time().timeName(),
mesh()
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh(),
dimensionedScalar("0", dimVolume/dimTime, 0.0)

View file

@ -412,7 +412,10 @@ tmp<surfaceScalarField> steadyStateDdtScheme<Type>::meshPhi
(
"meshPhi",
mesh().time().timeName(),
mesh()
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh(),
dimensionedScalar("0", dimVolume/dimTime, 0.0)

View file

@ -46,15 +46,6 @@ wedgeFvPatch::wedgeFvPatch(const polyPatch& patch, const fvBoundaryMesh& bm)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::vectorField> Foam::wedgeFvPatch::delta() const
{
const vectorField nHat(nf());
return nHat*(nHat & (Cf() - Cn()));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View file

@ -83,9 +83,6 @@ public:
{
return wedgePolyPatch_.cellT();
}
//- Return cell-centre to face normal vector
virtual tmp<vectorField> delta() const;
};

View file

@ -140,7 +140,9 @@ const Foam::scalarField& Foam::fvPatch::magSf() const
Foam::tmp<Foam::vectorField> Foam::fvPatch::delta() const
{
return Cf() - Cn();
// Use patch-normal delta for all non-coupled BCs
const vectorField nHat(nf());
return nHat*(nHat & (Cf() - Cn()));
}

View file

@ -103,8 +103,12 @@ bool Foam::trackedParticle::move
scalar tEnd = (1.0 - stepFraction())*trackTime;
scalar dtMax = tEnd;
if (tEnd <= SMALL)
if (tEnd <= SMALL && onBoundary())
{
// This is a hack to handle particles reaching their endpoint
// on a processor boundary. If the endpoint is on a processor face
// it currently gets transferred backwards and forwards infinitely.
// Remove the particle
td.keepParticle = false;
}

View file

@ -30,7 +30,7 @@ License
template<class SourcePatch, class TargetPatch>
void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
(
boolList& mapFlag,
labelList& mapFlag,
labelList& srcTgtSeed,
DynamicList<label>& srcSeeds,
DynamicList<label>& nonOverlapFaces,
@ -50,19 +50,32 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
{
label srcI = srcNbr[i];
if (mapFlag[srcI] && srcTgtSeed[srcI] == -1)
if ((mapFlag[srcI] == 0) && (srcTgtSeed[srcI] == -1))
{
// first attempt: match by comparing face centres
const face& srcF = this->srcPatch_[srcI];
const vector srcN = srcF.normal(srcPoints);
const point& srcC = srcCf[srcI];
scalar tol = GREAT;
forAll(srcF, fpI)
{
const point& p = srcPoints[srcF[fpI]];
scalar d2 = magSqr(p - srcC);
if (d2 < tol)
{
tol = d2;
}
}
tol = max(SMALL, 0.0001*sqrt(tol));
bool found = false;
forAll(tgtNbr, j)
{
label tgtI = tgtNbr[j];
const face& tgtF = this->tgtPatch_[tgtI];
pointHit ray = tgtF.ray(srcCf[srcI], srcN, tgtPoints);
const point tgtC = tgtF.centre(tgtPoints);
if (ray.hit())
if (mag(srcC - tgtC) < tol)
{
// new match - append to lists
found = true;
@ -74,11 +87,57 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
}
}
// second attempt: match by shooting a ray into the tgt face
if (!found)
{
// no match available for source face srcI
mapFlag[srcI] = false;
const vector srcN = srcF.normal(srcPoints);
forAll(tgtNbr, j)
{
label tgtI = tgtNbr[j];
const face& tgtF = this->tgtPatch_[tgtI];
pointHit ray = tgtF.ray(srcCf[srcI], srcN, tgtPoints);
if (ray.hit())
{
// new match - append to lists
found = true;
srcTgtSeed[srcI] = tgtI;
srcSeeds.append(srcI);
break;
}
}
}
// no match available for source face srcI
if (!found)
{
mapFlag[srcI] = -1;
nonOverlapFaces.append(srcI);
if (debug)
{
Pout<< "source face not found: id=" << srcI
<< " centre=" << srcCf[srcI]
<< " normal=" << srcF.normal(srcPoints)
<< " points=" << srcF.points(srcPoints)
<< endl;
Pout<< "target neighbours:" << nl;
forAll(tgtNbr, j)
{
label tgtI = tgtNbr[j];
const face& tgtF = this->tgtPatch_[tgtI];
Pout<< "face id: " << tgtI
<< " centre=" << tgtF.centre(tgtPoints)
<< " normal=" << tgtF.normal(tgtPoints)
<< " points=" << tgtF.points(tgtPoints)
<< endl;
}
}
}
}
}
@ -96,6 +155,36 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
}
template<class SourcePatch, class TargetPatch>
void Foam::directAMI<SourcePatch, TargetPatch>::restartAdvancingFront
(
labelList& mapFlag,
DynamicList<label>& nonOverlapFaces,
label& srcFaceI,
label& tgtFaceI
) const
{
forAll(mapFlag, faceI)
{
if (mapFlag[faceI] == 0)
{
tgtFaceI = this->findTargetFace(faceI);
if (tgtFaceI < 0)
{
mapFlag[faceI] = -1;
nonOverlapFaces.append(faceI);
}
else
{
srcFaceI = faceI;
break;
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
@ -176,15 +265,19 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
srcTgtSeed[srcFaceI] = tgtFaceI;
// list to keep track of whether src face can be mapped
boolList mapFlag(srcAddr.size(), true);
// 1 = mapped, 0 = untested, -1 = cannot map
labelList mapFlag(srcAddr.size(), 0);
label nTested = 0;
DynamicList<label> nonOverlapFaces;
do
{
srcAddr[srcFaceI].append(tgtFaceI);
tgtAddr[tgtFaceI].append(srcFaceI);
mapFlag[srcFaceI] = false;
mapFlag[srcFaceI] = 1;
nTested++;
// Do advancing front starting from srcFaceI, tgtFaceI
appendToDirectSeeds
@ -196,6 +289,12 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
srcFaceI,
tgtFaceI
);
if (srcFaceI < 0 && nTested < this->srcPatch_.size())
{
restartAdvancingFront(mapFlag, nonOverlapFaces, srcFaceI, tgtFaceI);
}
} while (srcFaceI >= 0);
if (nonOverlapFaces.size() != 0)
@ -211,14 +310,16 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
forAll(srcAddr, i)
{
scalar magSf = this->srcMagSf_[i];
srcAddress[i].transfer(srcAddr[i]);
// srcWeights[i] = scalarList(srcAddr[i].size(), magSf);
srcWeights[i] = scalarList(1, magSf);
srcAddress[i].transfer(srcAddr[i]);
}
forAll(tgtAddr, i)
{
scalar magSf = this->tgtMagSf_[i];
tgtAddress[i].transfer(tgtAddr[i]);
// tgtWeights[i] = scalarList(tgtAddr[i].size(), magSf);
tgtWeights[i] = scalarList(1, magSf);
tgtAddress[i].transfer(tgtAddr[i]);
}
}

View file

@ -67,7 +67,7 @@ private:
//- Append to list of src face seed indices
void appendToDirectSeeds
(
boolList& mapFlag,
labelList& mapFlag,
labelList& srcTgtSeed,
DynamicList<label>& srcSeeds,
DynamicList<label>& nonOverlapFaces,
@ -75,6 +75,16 @@ private:
label& tgtFaceI
) const;
//- Restart the advancing front - typically happens for
// disconnected regions
void restartAdvancingFront
(
labelList& mapFlag,
DynamicList<label>& nonOverlapFaces,
label& srcFaceI,
label& tgtFaceI
) const;
// Evaluation

View file

@ -99,11 +99,10 @@ void Foam::fieldValues::faceSource::setFaceZoneFaces()
const faceZone& fZone = mesh().faceZones()[zoneId];
faceId_.setSize(fZone.size());
facePatchId_.setSize(fZone.size());
faceSign_.setSize(fZone.size());
DynamicList<label> faceIds(fZone.size());
DynamicList<label> facePatchIds(fZone.size());
DynamicList<label> faceSigns(fZone.size());
label count = 0;
forAll(fZone, i)
{
label faceI = fZone[i];
@ -145,27 +144,26 @@ void Foam::fieldValues::faceSource::setFaceZoneFaces()
{
if (fZone.flipMap()[i])
{
faceSign_[count] = -1;
faceSigns.append(-1);
}
else
{
faceSign_[count] = 1;
faceSigns.append(1);
}
faceId_[count] = faceId;
facePatchId_[count] = facePatchId;
count++;
faceIds.append(faceId);
facePatchIds.append(facePatchId);
}
}
faceId_.setSize(count);
facePatchId_.setSize(count);
faceSign_.setSize(count);
faceId_.transfer(faceIds);
facePatchId_.transfer(facePatchIds);
faceSign_.transfer(faceSigns);
nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
if (debug)
{
Pout<< "Original face zone size = " << fZone.size()
<< ", new size = " << count << endl;
<< ", new size = " << faceId_.size() << endl;
}
}
@ -445,12 +443,43 @@ void Foam::fieldValues::faceSource::initialise(const dictionary& dict)
if (dict.readIfPresent("weightField", weightFieldName_))
{
Info<< " weight field = " << weightFieldName_;
Info<< " weight field = " << weightFieldName_ << nl;
}
if (dict.found("orientedWeightField"))
{
if (weightFieldName_ == "none")
{
dict.lookup("orientedWeightField") >> weightFieldName_;
Info<< " weight field = " << weightFieldName_ << nl;
orientWeightField_ = true;
}
else
{
FatalIOErrorIn
(
"void Foam::fieldValues::faceSource::initialise"
"("
"const dictionary&"
")",
dict
)
<< "Either weightField or orientedWeightField can be supplied, "
<< "but not both"
<< exit(FatalIOError);
}
}
List<word> orientedFields;
if (dict.readIfPresent("orientedFields", orientedFields))
{
orientedFieldsStart_ = fields_.size();
fields_.append(orientedFields);
}
if (dict.readIfPresent("scaleFactor", scaleFactor_))
{
Info<< " scale factor = " << scaleFactor_;
Info<< " scale factor = " << scaleFactor_ << nl;
}
Info<< nl << endl;
@ -581,6 +610,8 @@ Foam::fieldValues::faceSource::faceSource
source_(sourceTypeNames_.read(dict.lookup("source"))),
operation_(operationTypeNames_.read(dict.lookup("operation"))),
weightFieldName_("none"),
orientWeightField_(false),
orientedFieldsStart_(labelMax),
scaleFactor_(1.0),
nFaces_(0),
faceId_(),
@ -628,24 +659,42 @@ void Foam::fieldValues::faceSource::write()
totalArea = gSum(filterField(mesh().magSf(), false));
}
if (Pstream::master())
{
file() << obr_.time().value() << tab << totalArea;
}
// construct weight field
scalarField weightField(faceId_.size(), 1.0);
if (weightFieldName_ != "none")
{
weightField =
getFieldValues<scalar>
(
weightFieldName_,
true,
orientWeightField_
);
}
// Combine onto master
combineFields(weightField);
// process the fields
forAll(fields_, i)
{
const word& fieldName = fields_[i];
bool processed = false;
bool ok = false;
processed = processed || writeValues<scalar>(fieldName);
processed = processed || writeValues<vector>(fieldName);
processed = processed || writeValues<sphericalTensor>(fieldName);
processed = processed || writeValues<symmTensor>(fieldName);
processed = processed || writeValues<tensor>(fieldName);
bool orient = i >= orientedFieldsStart_;
ok = ok || writeValues<scalar>(fieldName, weightField, orient);
ok = ok || writeValues<vector>(fieldName, weightField, orient);
ok = ok
|| writeValues<sphericalTensor>(fieldName, weightField, orient);
ok = ok || writeValues<symmTensor>(fieldName, weightField, orient);
ok = ok || writeValues<tensor>(fieldName, weightField, orient);
if (!processed)
if (!ok)
{
WarningIn("void Foam::fieldValues::faceSource::write()")
<< "Requested field " << fieldName

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -71,8 +71,10 @@ Description
sourceName | name of face source if required | no |
operation | operation to perform | yes |
weightField | name of field to apply weighting | no |
orientedWeightField | name of oriented field to apply weighting | no |
scaleFactor | scale factor | no | 1
fields | list of fields to operate on | yes |
orientedFields | list of oriented fields to operate on | no |
\endtable
\linebreak
@ -109,8 +111,8 @@ Note
- faces on empty patches get ignored
- if the field is a volField the \c faceZone can only consist of boundary
faces
- all fields are oriented according to the \c faceZone (so you might
e.g. see negative pressure)
- the `oriented' entries relate to mesh-oriented fields, such as the
flux, phi. These fields will be oriented according to the face normals.
- using \c sampledSurfaces:
- not available for surface fields
- if interpolate=true they use \c interpolationCellPoint
@ -128,6 +130,7 @@ SeeAlso
SourceFiles
faceSource.C
faceSourceTemplates.C
\*---------------------------------------------------------------------------*/
@ -242,6 +245,12 @@ protected:
//- Weight field name - optional
word weightFieldName_;
//- Flag to indicate if flipMap should be applied to the weight field
bool orientWeightField_;
//- Start index of fields that require application of flipMap
label orientedFieldsStart_;
//- Scale factor - optional
scalar scaleFactor_;
@ -282,7 +291,8 @@ protected:
tmp<Field<Type> > getFieldValues
(
const word& fieldName,
const bool mustGet = false
const bool mustGet = false,
const bool applyOrientation = false
) const;
//- Apply the 'operation' to the values. Operation has to
@ -356,7 +366,12 @@ public:
//- Templated helper function to output field values
template<class Type>
bool writeValues(const word& fieldName);
bool writeValues
(
const word& fieldName,
const scalarField& weightField,
const bool orient
);
//- Filter a surface field according to faceIds
template<class Type>

View file

@ -54,7 +54,8 @@ template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::getFieldValues
(
const word& fieldName,
const bool mustGet
const bool mustGet,
const bool applyOrientation
) const
{
typedef GeometricField<Type, fvsPatchField, surfaceMesh> sf;
@ -62,7 +63,7 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::getFieldValues
if (source_ != stSampledSurface && obr_.foundObject<sf>(fieldName))
{
return filterField(obr_.lookupObject<sf>(fieldName), true);
return filterField(obr_.lookupObject<sf>(fieldName), applyOrientation);
}
else if (obr_.foundObject<vf>(fieldName))
{
@ -103,7 +104,7 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::getFieldValues
}
else
{
return filterField(fld, true);
return filterField(fld, applyOrientation);
}
}
@ -115,6 +116,7 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::getFieldValues
"Foam::fieldValues::faceSource::getFieldValues"
"("
"const word&, "
"const bool, "
"const bool"
") const"
) << "Field " << fieldName << " not found in database"
@ -267,22 +269,20 @@ Type Foam::fieldValues::faceSource::processValues
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
bool Foam::fieldValues::faceSource::writeValues(const word& fieldName)
bool Foam::fieldValues::faceSource::writeValues
(
const word& fieldName,
const scalarField& weightField,
const bool orient
)
{
const bool ok = validField<Type>(fieldName);
if (ok)
{
Field<Type> values(getFieldValues<Type>(fieldName));
scalarField weightField(values.size(), 1.0);
if (weightFieldName_ != "none")
{
weightField = getFieldValues<scalar>(weightFieldName_, true);
}
Field<Type> values(getFieldValues<Type>(fieldName, true, orient));
vectorField Sf;
if (surfacePtr_.valid())
{
// Get oriented Sf
@ -291,13 +291,12 @@ bool Foam::fieldValues::faceSource::writeValues(const word& fieldName)
else
{
// Get oriented Sf
Sf = filterField(mesh().Sf(), false);
Sf = filterField(mesh().Sf(), true);
}
// Combine onto master
combineFields(values);
combineFields(Sf);
combineFields(weightField);
// Write raw values on surface if specified
if (surfaceWriterPtr_.valid())
@ -378,7 +377,8 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::filterField
(
"fieldValues::faceSource::filterField"
"("
"const GeometricField<Type, fvPatchField, volMesh>&"
"const GeometricField<Type, fvPatchField, volMesh>&, "
"const bool"
") const"
) << type() << " " << name_ << ": "
<< sourceTypeNames_[source_] << "(" << sourceName_ << "):"

View file

@ -52,7 +52,7 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
os.writeKeyword("centreOfMass")
<< initialCentreOfMass_ << token::END_STATEMENT << nl;
os.writeKeyword("orientation")
os.writeKeyword("initialOrientation")
<< initialQ_ << token::END_STATEMENT << nl;
os.writeKeyword("mass")
<< mass_ << token::END_STATEMENT << nl;

View file

@ -33,6 +33,8 @@ Foam::PengRobinsonGas<Specie>::PengRobinsonGas(Istream& is)
:
Specie(is),
Tc_(readScalar(is)),
Vc_(readScalar(is)),
Zc_(readScalar(is)),
Pc_(readScalar(is)),
omega_(readScalar(is))
{
@ -48,9 +50,13 @@ Foam::PengRobinsonGas<Specie>::PengRobinsonGas
:
Specie(dict),
Tc_(readScalar(dict.subDict("equationOfState").lookup("Tc"))),
Vc_(readScalar(dict.subDict("equationOfState").lookup("Vc"))),
Zc_(1.0),
Pc_(readScalar(dict.subDict("equationOfState").lookup("Pc"))),
omega_(readScalar(dict.subDict("equationOfState").lookup("omega")))
{}
{
Zc_ = Pc_*Vc_/(specie::RR*Tc_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -74,6 +80,8 @@ Foam::Ostream& Foam::operator<<
{
os << static_cast<const Specie&>(pg)
<< token::SPACE << pg.Tc_
<< token::SPACE << pg.Vc_
<< token::SPACE << pg.Zc_
<< token::SPACE << pg.Pc_
<< token::SPACE << pg.omega_;

View file

@ -98,10 +98,16 @@ class PengRobinsonGas
//- Critical Temperature [K]
scalar Tc_;
//- Critical volume [m^3/kmol]
scalar Vc_;
//- Critical compression factor [-]
scalar Zc_;
//- Critical Pressure [Pa]
scalar Pc_;
//- Accentric factor [-]
//- Acentric factor [-]
scalar omega_;
@ -114,6 +120,8 @@ public:
(
const Specie& sp,
const scalar& Tc,
const scalar& Vc,
const scalar& Zc,
const scalar& Pc,
const scalar& omega
);
@ -163,7 +171,7 @@ public:
//- Return compressibility rho/p [s^2/m^2]
inline scalar psi(scalar p, scalar T) const;
//- Return compression factor []
//- Return compression factor [-]
inline scalar Z(scalar p, scalar T) const;
//- Return (cp - cv) [J/(kmol K]

View file

@ -21,7 +21,6 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "PengRobinsonGas.H"
@ -34,12 +33,16 @@ inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
(
const Specie& sp,
const scalar& Tc,
const scalar& Vc,
const scalar& Zc,
const scalar& Pc,
const scalar& omega
)
:
Specie(sp),
Tc_(Tc),
Vc_(Vc),
Zc_(Zc),
Pc_(Pc),
omega_(omega)
{}
@ -55,9 +58,11 @@ inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
)
:
Specie(name, pg),
Tc_(pg.Tc),
Pc_(pg.Pc),
omega_(pg.omega)
Tc_(pg.Tc_),
Pc_(pg.Pc_),
Vc_(pg.Vc_),
Zc_(pg.Zc_),
omega_(pg.omega_)
{}
@ -96,6 +101,7 @@ Foam::PengRobinsonGas<Specie>::New
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Specie>
@ -214,7 +220,9 @@ inline void Foam::PengRobinsonGas<Specie>::operator+=
scalar molr2 = pg.nMoles()/this->nMoles();
Tc_ = molr1*Tc_ + molr2*pg.Tc_;
Pc_ = molr1*Pc_ + molr2*pg.Pc_;
Vc_ = molr1*Vc_ + molr2*pg.Vc_;
Zc_ = molr1*Zc_ + molr2*pg.Zc_;
Pc_ = specie::RR*Zc_*Tc_/Vc_;
omega_ = molr1*omega_ + molr2*pg.omega_;
}
@ -233,7 +241,9 @@ inline void Foam::PengRobinsonGas<Specie>::operator-=
scalar molr2 = pg.nMoles()/this->nMoles();
Tc_ = molr1*Tc_ - molr2*pg.Tc_;
Pc_ = molr1*Pc_ - molr2*pg.Pc_;
Vc_ = molr1*Vc_ - molr2*pg.Vc_;
Zc_ = molr1*Zc_ - molr2*pg.Zc_;
Pc_ = specie::RR*Zc_*Tc_/Vc_;
omega_ = molr1*omega_ - molr2*pg.omega_;
}
@ -259,12 +269,18 @@ Foam::PengRobinsonGas<Specie> Foam::operator+
scalar molr1 = pg1.nMoles()/nMoles;
scalar molr2 = pg2.nMoles()/nMoles;
const scalar Tc = molr1*pg1.Tc_ + molr2*pg2.Tc_;
const scalar Vc = molr1*pg1.Vc_ + molr2*pg2.Vc_;
const scalar Zc = molr1*pg1.Zc_ + molr2*pg2.Zc_;
return PengRobinsonGas<Specie>
(
static_cast<const Specie&>(pg1)
+ static_cast<const Specie&>(pg2),
molr1*pg1.Tc_ + molr2*pg2.Tc_,
molr1*pg1.Pc_ + molr2*pg2.Pc_,
Tc,
Vc,
Zc,
specie::RR*Zc*Tc/Vc,
molr1*pg1.omega_ + molr2*pg2.omega_
);
}
@ -281,12 +297,18 @@ Foam::PengRobinsonGas<Specie> Foam::operator-
scalar molr1 = pg1.nMoles()/nMoles;
scalar molr2 = pg2.nMoles()/nMoles;
const scalar Tc = molr1*pg1.Tc_ + molr2*pg2.Tc_;
const scalar Vc = molr1*pg1.Vc_ + molr2*pg2.Vc_;
const scalar Zc = molr1*pg1.Zc_ + molr2*pg2.Zc_;
return PengRobinsonGas<Specie>
(
static_cast<const Specie&>(pg1)
- static_cast<const Specie&>(pg2),
molr1*pg1.Tc_ - molr2*pg2.Tc_,
molr1*pg1.Pc_ - molr2*pg2.Pc_,
Tc,
Vc,
Zc,
specie::RR*Zc*Tc/Vc,
molr1*pg1.omega_ - molr2*pg2.omega_
);
}
@ -303,6 +325,8 @@ Foam::PengRobinsonGas<Specie> Foam::operator*
(
s*static_cast<const Specie&>(pg),
pg.Tc_,
pg.Vc_,
pg.Zc_,
pg.Pc_,
pg.omega_
);