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:
Henry
2011-04-20 18:30:48 +01:00
parent 119f1fdc62 ba01d338fb
commit 188469b994
23 changed files with 345 additions and 315 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
View File

@ -1210,12 +1210,13 @@ int main(int argc, char *argv[])
<< endl;
label nSide = 0;
forAll(zoneSidePatch, zoneI)
{
if (oneD)
{
// Always add empty patches, one per zone.
word patchName = faceZones[zoneI].name() + "_" + "side";
// Reuse single empty patch.
word patchName = "oneDEmptPatch";
zoneSidePatch[zoneI] = addPatch<emptyPolyPatch>
(

58
applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
View File

@ -6,6 +6,7 @@
#include "faceSet.H"
#include "pointSet.H"
#include "IOmanip.H"
#include "emptyPolyPatch.H"
Foam::label Foam::checkTopology
(
@ -21,6 +22,29 @@ Foam::label Foam::checkTopology
// Check if the boundary definition is unique
mesh.boundaryMesh().checkDefinition(true);
// Check that empty patches cover all sides of the mesh
{
label nEmpty = 0;
forAll(mesh.boundaryMesh(), patchI)
{
if (isA<emptyPolyPatch>(mesh.boundaryMesh()[patchI]))
{
nEmpty += mesh.boundaryMesh()[patchI].size();
}
}
reduce(nEmpty, sumOp<label>());
label nTotCells = returnReduce(mesh.cells().size(), sumOp<label>());
// These are actually warnings, not errors.
if (nEmpty % nTotCells)
{
Info<< " ***Total number of faces on empty patches"
<< " is not divisible by the number of cells in the mesh."
<< " Hence this mesh is not 1D or 2D."
<< endl;
}
}
// Check if the boundary processor patches are correct
mesh.boundaryMesh().checkParallelSync(true);
@ -41,6 +65,8 @@ Foam::label Foam::checkTopology
noFailedChecks++;
}
{
pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
if (mesh.checkPoints(true, &points))
@ -74,6 +100,22 @@ Foam::label Foam::checkTopology
}
}
{
faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
if (mesh.checkFaceVertices(true, &faces))
{
noFailedChecks++;
label nFaces = returnReduce(faces.size(), sumOp<label>());
Info<< " <<Writing " << nFaces
<< " faces with out-of-range or duplicate vertices to set "
<< faces.name() << endl;
faces.instance() = mesh.pointsInstance();
faces.write();
}
}
if (allTopology)
{
cellSet cells(mesh, "zipUpCells", mesh.nCells()/100);
@ -91,22 +133,6 @@ Foam::label Foam::checkTopology
}
}
{
faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
if (mesh.checkFaceVertices(true, &faces))
{
noFailedChecks++;
label nFaces = returnReduce(faces.size(), sumOp<label>());
Info<< " <<Writing " << nFaces
<< " faces with out-of-range or duplicate vertices to set "
<< faces.name() << endl;
faces.instance() = mesh.pointsInstance();
faces.write();
}
}
if (allTopology)
{
faceSet faces(mesh, "edgeFaces", mesh.nFaces()/100);

18
src/OpenFOAM/db/Time/Time.C
View File

@ -228,6 +228,8 @@ Foam::Time::Time
objectRegistry(*this),
libs_(),
controlDict_
(
IOobject
@ -257,9 +259,10 @@ Foam::Time::Time
graphFormat_("raw"),
runTimeModifiable_(true),
libs_(controlDict_, "libs"),
functionObjects_(*this)
{
libs_.open(controlDict_, "libs");
// Explicitly set read flags on objectRegistry so anything constructed
// from it reads as well (e.g. fvSolution).
readOpt() = IOobject::MUST_READ_IF_MODIFIED;
@ -313,6 +316,8 @@ Foam::Time::Time
objectRegistry(*this),
libs_(),
controlDict_
(
IOobject
@ -343,9 +348,11 @@ Foam::Time::Time
graphFormat_("raw"),
runTimeModifiable_(true),
libs_(controlDict_, "libs"),
functionObjects_(*this)
{
libs_.open(controlDict_, "libs");
// Explicitly set read flags on objectRegistry so anything constructed
// from it reads as well (e.g. fvSolution).
readOpt() = IOobject::MUST_READ_IF_MODIFIED;
@ -401,6 +408,8 @@ Foam::Time::Time
objectRegistry(*this),
libs_(),
controlDict_
(
IOobject
@ -430,9 +439,10 @@ Foam::Time::Time
graphFormat_("raw"),
runTimeModifiable_(true),
libs_(controlDict_, "libs"),
functionObjects_(*this)
{}
{
libs_.open(controlDict_, "libs");
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //

7
src/OpenFOAM/db/Time/Time.H
View File

@ -75,6 +75,10 @@ class Time
//- file-change monitor for all registered files
mutable autoPtr<fileMonitor> monitorPtr_;
//- Any loaded dynamic libraries. Make sure to construct before
// reading controlDict.
dlLibraryTable libs_;
//- The controlDict
IOdictionary controlDict_;
@ -166,9 +170,6 @@ private:
//- Is runtime modification of dictionaries allowed?
Switch runTimeModifiable_;
//- Any loaded dynamic libraries
dlLibraryTable libs_;
//- Function objects executed at start and on ++, +=
mutable functionObjectList functionObjects_;

18
src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C
View File

@ -131,11 +131,7 @@ bool Foam::functionEntries::codeStream::execute
// see if library is loaded
void* lib = NULL;
if
(
isA<IOdictionary>(topDict(parentDict))
&& parentDict.dictName() != Time::controlDictName
)
if (isA<IOdictionary>(topDict(parentDict)))
{
lib = libs(parentDict).findLibrary(libPath);
}
@ -150,11 +146,7 @@ bool Foam::functionEntries::codeStream::execute
// avoid compilation if possible by loading an existing library
if (!lib)
{
if
(
isA<IOdictionary>(topDict(parentDict))
&& parentDict.dictName() != Time::controlDictName
)
if (isA<IOdictionary>(topDict(parentDict)))
{
// Cached access to dl libs. Guarantees clean up upon destruction
// of Time.
@ -223,11 +215,7 @@ bool Foam::functionEntries::codeStream::execute
// all processes must wait for compile to finish
reduce(create, orOp<bool>());
if
(
isA<IOdictionary>(topDict(parentDict))
&& parentDict.dictName() != Time::controlDictName
)
if (isA<IOdictionary>(topDict(parentDict)))
{
// Cached access to dl libs. Guarantees clean up upon destruction
// of Time.

24
src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C
View File

@ -25,6 +25,12 @@ License
#include "dlLibraryTable.H"
#include "OSspecific.H"
#include "long.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::dlLibraryTable, 0);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -55,7 +61,11 @@ Foam::dlLibraryTable::~dlLibraryTable()
// bug in dlclose - does not call static destructors of
// loaded library when actually unloading the library.
// See https://bugzilla.novell.com/show_bug.cgi?id=680125 and 657627.
// Seems related to using a non-system compiler!
if (debug)
{
Info<< "dlLibraryTable::~dlLibraryTable() : closing " << iter()
<< " with handle " << long(iter.key()) << endl;
}
dlClose(iter.key());
}
}
@ -73,6 +83,12 @@ bool Foam::dlLibraryTable::open
{
void* functionLibPtr = dlOpen(functionLibName);
if (debug)
{
Info<< "dlLibraryTable::open : opened " << functionLibName
<< " resulting in handle " << long(functionLibPtr) << endl;
}
if (!functionLibPtr)
{
if (verbose)
@ -107,6 +123,12 @@ bool Foam::dlLibraryTable::close
void* libPtr = findLibrary(functionLibName);
if (libPtr)
{
if (debug)
{
Info<< "dlLibraryTable::close : closing " << functionLibName
<< " with handle " << long(libPtr) << endl;
}
erase(libPtr);
if (!dlClose(libPtr))

3
src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.H
View File

@ -63,6 +63,9 @@ class dlLibraryTable
public:
// Declare name of the class and its debug switch
ClassName("dlLibraryTable");
// Constructors
//- Construct null

73
src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
View File

@ -81,30 +81,6 @@ void Foam::addPatchCellLayer::addVertex
f[fp++] = pointI;
}
}
// Check for duplicates.
if (debug)
{
label n = 0;
for (label i = 0; i < fp; i++)
{
if (f[i] == pointI)
{
n++;
if (n == 2)
{
f.setSize(fp);
FatalErrorIn
(
"addPatchCellLayer::addVertex(const label, face&"
", label&)"
) << "Point " << pointI << " already present in face "
<< f << abort(FatalError);
}
}
}
}
}
@ -1641,11 +1617,11 @@ void Foam::addPatchCellLayer::setRefinement
{
label offset =
addedPoints_[vStart].size() - numEdgeSideFaces;
for (label ioff = offset; ioff > 0; ioff--)
for (label ioff = offset-1; ioff >= 0; ioff--)
{
addVertex
(
addedPoints_[vStart][ioff-1],
addedPoints_[vStart][ioff],
newFace,
newFp
);
@ -1660,6 +1636,51 @@ void Foam::addPatchCellLayer::setRefinement
newFace.setSize(newFp);
if (debug)
{
labelHashSet verts(2*newFace.size());
forAll(newFace, fp)
{
if (!verts.insert(newFace[fp]))
{
FatalErrorIn
(
"addPatchCellLayer::setRefinement(..)"
) << "Duplicate vertex in face"
<< " to be added." << nl
<< "newFace:" << newFace << nl
<< "points:"
<< UIndirectList<point>
(
meshMod.points(),
newFace
) << nl
<< "Layer:" << i
<< " out of:" << numEdgeSideFaces << nl
<< "ExtrudeEdge:" << meshEdgeI
<< " at:"
<< mesh_.edges()[meshEdgeI].line
(
mesh_.points()
) << nl
<< "string:" << stringedVerts
<< "stringpoints:"
<< UIndirectList<point>
(
pp.localPoints(),
stringedVerts
) << nl
<< "stringNLayers:"
<< UIndirectList<label>
(
nPointLayers,
stringedVerts
) << nl
<< abort(FatalError);
}
}
}
label nbrFaceI = nbrFace
(
pp.edgeFaces(),

128
src/finiteVolume/cfdTools/general/pimpleLoop/pimpleLoop.H
View File

@ -1,128 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::pimpleLoop
Description
PIMPLE loop class to formalise the iteration and automate the handling
of the "finalIteration" mesh data entry.
\*---------------------------------------------------------------------------*/
#ifndef pimpleLoop_H
#define pimpleLoop_H
#include "fvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class pimpleLoop Declaration
\*---------------------------------------------------------------------------*/
class pimpleLoop
{
// Private data
//- Reference to the mesh
fvMesh& mesh_;
//- Number of PIMPLE correctors
const int nCorr_;
//- Current PIMPLE corrector
int corr_;
// Private Member Functions
//- Disallow default bitwise copy construct
pimpleLoop(const pimpleLoop&);
//- Disallow default bitwise assignment
void operator=(const pimpleLoop&);
public:
// Constructors
//- Construct from components
pimpleLoop(fvMesh& mesh, const int nCorr)
:
mesh_(mesh),
nCorr_(nCorr),
corr_(0)
{}
//- Destructor
~pimpleLoop()
{}
// Member Functions
bool loop()
{
if (finalIter())
{
mesh_.data::add("finalIteration", true);
}
return corr_ < nCorr_;
}
bool finalIter() const
{
return corr_ == nCorr_-1;
}
// Member Operators
void operator++(int)
{
if (finalIter())
{
mesh_.data::remove("finalIteration");
}
corr_++;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
View File

@ -137,7 +137,7 @@ Foam::pimpleControl::pimpleControl(fvMesh& mesh)
{
Info<< " field " << residualControl_[i].name << token::TAB
<< ": relTol " << residualControl_[i].relTol
<< ", absTol " << residualControl_[i].absTol
<< ", tolerance " << residualControl_[i].absTol
<< nl;
}
Info<< endl;

4
src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C
View File

@ -68,7 +68,7 @@ bool Foam::simpleControl::criteriaSatisfied()
{
Info<< algorithmName_ << " solution statistics:" << endl;
Info<< " " << variableName << ": abs tol = " << residual
Info<< " " << variableName << ": tolerance = " << residual
<< " (" << residualControl_[fieldI].absTol << ")"
<< endl;
}
@ -96,7 +96,7 @@ Foam::simpleControl::simpleControl(fvMesh& mesh)
forAll(residualControl_, i)
{
Info<< " field " << residualControl_[i].name << token::TAB
<< " absTol " << residualControl_[i].absTol
<< " tolerance " << residualControl_[i].absTol
<< nl;
}
Info<< endl;

2
src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
View File

@ -68,7 +68,7 @@ void Foam::solutionControl::read(const bool absTolOnly)
if (iter().isDict())
{
const dictionary& fieldDict(iter().dict());
fd.absTol = readScalar(fieldDict.lookup("absTol"));
fd.absTol = readScalar(fieldDict.lookup("tolerance"));
fd.relTol = readScalar(fieldDict.lookup("relTol"));
fd.initialResidual = 0.0;
}

29
src/finiteVolume/fields/fvPatchFields/constraint/empty/emptyFvPatchField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -138,18 +138,21 @@ emptyFvPatchField<Type>::emptyFvPatchField
template<class Type>
void emptyFvPatchField<Type>::updateCoeffs()
{
if
(
this->patch().patch().size()
% this->dimensionedInternalField().mesh().nCells()
)
{
FatalErrorIn("emptyFvPatchField<Type>::updateCoeffs()")
<< "This mesh contains patches of type empty but is not 1D or 2D\n"
" by virtue of the fact that the number of faces of this\n"
" empty patch is not divisible by the number of cells."
<< exit(FatalError);
}
//- Check moved to checkMesh. Test here breaks down if multiple empty
// patches.
//if
//(
// this->patch().patch().size()
// % this->dimensionedInternalField().mesh().nCells()
//)
//{
// FatalErrorIn("emptyFvPatchField<Type>::updateCoeffs()")
// << "This mesh contains patches of type empty but is not"
// << "1D or 2D\n"
// " by virtue of the fact that the number of faces of this\n"
// " empty patch is not divisible by the number of cells."
// << exit(FatalError);
//}
}

2
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
View File

@ -518,7 +518,7 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::Re
const scalar muc
) const
{
return rhoc*mag(U - Uc_)*d/muc;
return rhoc*mag(U - Uc_)*d/(muc + ROOTVSMALL);
}

90
src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
View File

@ -63,7 +63,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::CpEff
template<class ParcelType>
template<class TrackData>
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HsEff
(
TrackData& td,
const scalar p,
@ -74,9 +74,9 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
) const
{
return
this->Y_[GAS]*td.cloud().composition().H(idG, YGas_, p, T)
+ this->Y_[LIQ]*td.cloud().composition().H(idL, YLiquid_, p, T)
+ this->Y_[SLD]*td.cloud().composition().H(idS, YSolid_, p, T);
this->Y_[GAS]*td.cloud().composition().Hs(idG, YGas_, p, T)
+ this->Y_[LIQ]*td.cloud().composition().Hs(idL, YLiquid_, p, T)
+ this->Y_[SLD]*td.cloud().composition().Hs(idS, YSolid_, p, T);
}
@ -326,7 +326,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
// Heat transfer
// ~~~~~~~~~~~~~
@ -368,7 +367,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
d0,
U0,
rho0,
0.5*(mass0 + mass1),
mass0,
Su,
dUTrans,
Spu
@ -383,25 +382,41 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Transfer mass lost from particle to carrier mass source
forAll(YGas_, i)
{
scalar dm = np0*dMassGas[i];
label gid = composition.localToGlobalCarrierId(GAS, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
scalar hs = composition.carrier().Hs(gid, T0);
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().UTrans()[cellI] += dm*U0;
td.cloud().hsTrans()[cellI] += dm*hs;
}
forAll(YLiquid_, i)
{
scalar dm = np0*dMassLiquid[i];
label gid = composition.localToGlobalCarrierId(LIQ, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
scalar hs = composition.carrier().Hs(gid, T0);
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().UTrans()[cellI] += dm*U0;
td.cloud().hsTrans()[cellI] += dm*hs;
}
/*
// No mapping between solid components and carrier phase
forAll(YSolid_, i)
{
scalar dm = np0*dMassSolid[i];
label gid = composition.localToGlobalCarrierId(SLD, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
scalar hs = composition.carrier().Hs(gid, T0);
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().UTrans()[cellI] += dm*U0;
td.cloud().hsTrans()[cellI] += dm*hs;
}
*/
forAll(dMassSRCarrier, i)
{
td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
scalar dm = np0*dMassSRCarrier[i];
scalar hs = composition.carrier().Hs(i, T0);
td.cloud().rhoTrans(i)[cellI] += dm;
td.cloud().UTrans()[cellI] += dm*U0;
td.cloud().hsTrans()[cellI] += dm*hs;
}
// Update momentum transfer
@ -421,36 +436,38 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Remove the particle when mass falls below minimum threshold
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (mass1 < td.cloud().constProps().minParticleMass())
if (np0*mass1 < td.cloud().constProps().minParticleMass())
{
td.keepParticle = false;
if (td.cloud().solution().coupled())
{
scalar dm = np0*mass1;
// Absorb parcel into carrier phase
forAll(YGas_, i)
{
label gid = composition.localToGlobalCarrierId(GAS, i);
td.cloud().rhoTrans(gid)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
td.cloud().rhoTrans(gid)[cellI] += dm*YMix[GAS]*YGas_[i];
}
forAll(YLiquid_, i)
{
label gid = composition.localToGlobalCarrierId(LIQ, i);
td.cloud().rhoTrans(gid)[cellI] +=
np0*mass1*YMix[LIQ]*YLiquid_[i];
td.cloud().rhoTrans(gid)[cellI] += dm*YMix[LIQ]*YLiquid_[i];
}
/*
// No mapping between solid components and carrier phase
forAll(YSolid_, i)
{
label gid = composition.localToGlobalCarrierId(SLD, i);
td.cloud().rhoTrans(gid)[cellI] +=
np0*mass1*YMix[SLD]*YSolid_[i];
td.cloud().rhoTrans(gid)[cellI] += dm*YMix[SLD]*YSolid_[i];
}
*/
td.cloud().UTrans()[cellI] += np0*mass1*U1;
td.cloud().UTrans()[cellI] += dm*U1;
// enthalpy transfer accounted for via change in mass fractions
td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T1, idG, idL, idS);
td.cloud().addToMassPhaseChange(dm);
}
}
@ -531,28 +548,35 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
Sh -= dMassTot*td.cloud().constProps().LDevol()/dt;
// Molar average molecular weight of carrier mix
const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
// Update molar emissions
forAll(dMassDV, i)
if (td.cloud().heatTransfer().BirdCorrection())
{
// Molar average molecular weight of carrier mix
const scalar Wc =
max(SMALL, this->rhoc_*specie::RR*this->Tc_/this->pc_);
// Note: hardcoded gaseous diffusivities for now
// TODO: add to carrier thermo
const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
const label id = composition.localToGlobalCarrierId(GAS, i);
const scalar Cp = composition.carrier().Cp(id, Ts);
const scalar W = composition.carrier().W(id);
const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
// Dab calc'd using API vapour mass diffusivity function
const scalar Dab =
3.6059e-3*(pow(1.8*Ts, 1.75))*sqrt(1.0/W + 1.0/Wc)/(this->pc_*beta);
forAll(dMassDV, i)
{
const label id = composition.localToGlobalCarrierId(GAS, i);
const scalar Cp = composition.carrier().Cp(id, Ts);
const scalar W = composition.carrier().W(id);
const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
N += Ni;
NCpW += Ni*Cp*W;
Cs[id] += Ni*d/(2.0*Dab);
}
// Dab calc'd using API vapour mass diffusivity function
const scalar Dab =
3.6059e-3*(pow(1.8*Ts, 1.75))
*sqrt(1.0/W + 1.0/Wc)
/(this->pc_*beta);
N += Ni;
NCpW += Ni*Cp*W;
Cs[id] += Ni*d/(2.0*Dab);
}
}
}

4
src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
View File

@ -134,9 +134,9 @@ private:
const label idS
) const;
//- Return the mixture effective enthalpy
//- Return the mixture effective sensible enthalpy
template<class TrackData>
scalar HEff
scalar HsEff
(
TrackData& td,
const scalar p,

84
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
View File

@ -213,9 +213,17 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
sumYiCbrtW += Ys[i]*cbrtW;
}
Cps = max(Cps, ROOTVSMALL);
rhos *= pc_/(specie::RR*T);
rhos = max(rhos, ROOTVSMALL);
mus /= sumYiSqrtW;
mus = max(mus, ROOTVSMALL);
kappas /= sumYiCbrtW;
kappas = max(kappas, ROOTVSMALL);
Prs = Cps*mus/kappas;
}
@ -335,7 +343,9 @@ void Foam::ReactingParcel<ParcelType>::calc
Res = this->Re(U0, d0, rhos, mus);
// Update particle component mass and mass fractions
scalar mass1 = updateMassFraction(mass0, dMassPC, Y_);
scalarField dMass(dMassPC);
scalar mass1 = updateMassFraction(mass0, dMass, Y_);
// Heat transfer
@ -379,7 +389,7 @@ void Foam::ReactingParcel<ParcelType>::calc
d0,
U0,
rho0,
0.5*(mass0 + mass1),
mass0,
Su,
dUTrans,
Spu
@ -390,11 +400,16 @@ void Foam::ReactingParcel<ParcelType>::calc
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (td.cloud().solution().coupled())
{
// Transfer mass lost from particle to carrier mass source
forAll(dMassPC, i)
// Transfer mass lost to carrier mass and enthalpy sources
forAll(dMass, i)
{
scalar dm = np0*dMass[i];
label gid = composition.localToGlobalCarrierId(0, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
scalar hs = composition.carrier().Hs(gid, T0);
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().UTrans()[cellI] += dm*U0;
td.cloud().hsTrans()[cellI] += dm*hs;
}
// Update momentum transfer
@ -413,21 +428,27 @@ void Foam::ReactingParcel<ParcelType>::calc
// Remove the particle when mass falls below minimum threshold
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (mass1 < td.cloud().constProps().minParticleMass())
if (np0*mass1 < td.cloud().constProps().minParticleMass())
{
td.keepParticle = false;
if (td.cloud().solution().coupled())
{
scalar dm = np0*mass1;
// Absorb parcel into carrier phase
forAll(Y_, i)
{
scalar dmi = dm*Y_[i];
label gid = composition.localToGlobalCarrierId(0, i);
td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
}
td.cloud().UTrans()[cellI] += np0*mass1*U1;
scalar hs = composition.carrier().Hs(gid, T1);
// enthalpy transfer accounted for via change in mass fractions
td.cloud().rhoTrans(gid)[cellI] += dmi;
td.cloud().hsTrans()[cellI] += dmi*hs;
}
td.cloud().UTrans()[cellI] += dm*U1;
td.cloud().addToMassPhaseChange(dm);
}
}
@ -514,33 +535,44 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
// Add to cumulative phase change mass
td.cloud().addToMassPhaseChange(this->nParticle_*dMassTot);
// Average molecular weight of carrier mix - assumes perfect gas
const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
forAll(YComponents, i)
forAll(dMassPC, i)
{
const label idc = composition.localToGlobalCarrierId(idPhase, i);
const label idl = composition.globalIds(idPhase)[i];
const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, T);
Sh -= dMassPC[i]*dh/dt;
}
// Update particle surface thermo properties
const scalar Dab =
composition.liquids().properties()[idl].D(pc_, Ts, Wc);
const scalar Cp = composition.carrier().Cp(idc, Ts);
const scalar W = composition.carrier().W(idc);
const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
// Update molar emissions
if (td.cloud().heatTransfer().BirdCorrection())
{
// Average molecular weight of carrier mix - assumes perfect gas
const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
// Molar flux of species coming from the particle (kmol/m^2/s)
N += Ni;
// Sum of Ni*Cpi*Wi of emission species
NCpW += Ni*Cp*W;
forAll(dMassPC, i)
{
const label idc = composition.localToGlobalCarrierId(idPhase, i);
const label idl = composition.globalIds(idPhase)[i];
// Concentrations of emission species
Cs[idc] += Ni*d/(2.0*Dab);
const scalar Cp = composition.carrier().Cp(idc, Ts);
const scalar W = composition.carrier().W(idc);
const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
const scalar Dab =
composition.liquids().properties()[idl].D(pc_, Ts, Wc);
// Molar flux of species coming from the particle (kmol/m^2/s)
N += Ni;
// Sum of Ni*Cpi*Wi of emission species
NCpW += Ni*Cp*W;
// Concentrations of emission species
Cs[idc] += Ni*d/(2.0*Dab);
}
}
}

3
src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
View File

@ -153,7 +153,10 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
mus = td.muInterp().interpolate(this->position(), tetIs)/TRatio;
Pr = td.cloud().constProps().Pr();
Pr = max(ROOTVSMALL, Pr);
kappas = Cpc_*mus/Pr;
kappas = max(ROOTVSMALL, kappas);
}

89
src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
View File

@ -129,7 +129,8 @@ Foam::refinementFeatures::refinementFeatures
}
Info<< "Detected " << featurePoints.size()
<< " featurePoints out of " << points.size() << endl;
<< " featurePoints out of " << points.size()
<< " on feature " << eMesh.name() << endl;
pointTrees_.set
(
@ -164,26 +165,30 @@ void Foam::refinementFeatures::findNearestEdge
forAll(edgeTrees_, featI)
{
const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
forAll(samples, sampleI)
if (tree.shapes().size() > 0)
{
const point& sample = samples[sampleI];
scalar distSqr;
if (nearInfo[sampleI].hit())
forAll(samples, sampleI)
{
distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
}
else
{
distSqr = nearestDistSqr[sampleI];
}
const point& sample = samples[sampleI];
pointIndexHit info = tree.findNearest(sample, distSqr);
scalar distSqr;
if (nearInfo[sampleI].hit())
{
distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
}
else
{
distSqr = nearestDistSqr[sampleI];
}
if (info.hit())
{
nearInfo[sampleI] = info;
nearFeature[sampleI] = featI;
pointIndexHit info = tree.findNearest(sample, distSqr);
if (info.hit())
{
nearInfo[sampleI] = info;
nearFeature[sampleI] = featI;
}
}
}
}
@ -206,33 +211,37 @@ void Foam::refinementFeatures::findNearestPoint
forAll(pointTrees_, featI)
{
const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
forAll(samples, sampleI)
if (tree.shapes().pointLabels().size() > 0)
{
const point& sample = samples[sampleI];
scalar distSqr;
if (nearFeature[sampleI] != -1)
forAll(samples, sampleI)
{
label nearFeatI = nearFeature[sampleI];
const indexedOctree<treeDataPoint>& nearTree =
pointTrees_[nearFeatI];
label featPointI =
nearTree.shapes().pointLabels()[nearIndex[sampleI]];
const point& featPt =
operator[](nearFeatI).points()[featPointI];
distSqr = magSqr(featPt-sample);
}
else
{
distSqr = nearestDistSqr[sampleI];
}
const point& sample = samples[sampleI];
pointIndexHit info = tree.findNearest(sample, distSqr);
scalar distSqr;
if (nearFeature[sampleI] != -1)
{
label nearFeatI = nearFeature[sampleI];
const indexedOctree<treeDataPoint>& nearTree =
pointTrees_[nearFeatI];
label featPointI =
nearTree.shapes().pointLabels()[nearIndex[sampleI]];
const point& featPt =
operator[](nearFeatI).points()[featPointI];
distSqr = magSqr(featPt-sample);
}
else
{
distSqr = nearestDistSqr[sampleI];
}
if (info.hit())
{
nearFeature[sampleI] = featI;
nearIndex[sampleI] = info.index();
pointIndexHit info = tree.findNearest(sample, distSqr);
if (info.hit())
{
nearFeature[sampleI] = featI;
nearIndex[sampleI] = info.index();
}
}
}
}

3
src/regionModels/regionModel/regionModel1D/regionModel1D.C
View File

@ -114,6 +114,7 @@ void Foam::regionModels::regionModel1D::initialise()
boundaryFaceCells_[localPyrolysisFaceI].transfer(cellIDs);
localPyrolysisFaceI++;
nLayers_ = nCells;
}
}
@ -266,6 +267,7 @@ Foam::regionModels::regionModel1D::regionModel1D(const fvMesh& mesh)
boundaryFaceFaces_(),
boundaryFaceCells_(),
boundaryFaceOppositeFace_(),
nLayers_(0),
nMagSfPtr_(NULL),
moveMesh_(false)
{}
@ -283,6 +285,7 @@ Foam::regionModels::regionModel1D::regionModel1D
boundaryFaceFaces_(regionMesh().nCells()),
boundaryFaceCells_(regionMesh().nCells()),
boundaryFaceOppositeFace_(regionMesh().nCells()),
nLayers_(0),
nMagSfPtr_(NULL),
moveMesh_(true)
{

6
src/regionModels/regionModel/regionModel1D/regionModel1D.H
View File

@ -88,6 +88,9 @@ protected:
//- Global boundary face IDs oppossite coupled patch
labelList boundaryFaceOppositeFace_;
//- Number of layers in the region
label nLayers_;
// Geometry
@ -157,6 +160,9 @@ public:
//- Return the face area magnitudes / [m2]
inline const surfaceScalarField& nMagSf() const;
//- Return the number of layers in the region
inline label nLayers() const;
};

6
src/regionModels/regionModel/regionModel1D/regionModel1DI.H
View File

@ -65,4 +65,10 @@ Foam::regionModels::regionModel1D::nMagSf() const
}
inline Foam::label Foam::regionModels::regionModel1D::nLayers() const
{
return nLayers_;
}
// ************************************************************************* //

2
tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSolution
View File

@ -63,7 +63,7 @@ PIMPLE
"(U|k|epsilon)"
{
relTol 0;
absTol 0.0001;
tolerance 0.0001;
}
}
}