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 /home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
mattijs
2013-09-03 14:42:34 +01:00
parent cadaee82b7 60f07d1f59
commit 2803571490
57 changed files with 798 additions and 556 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModels.C
View File

@ -1,44 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ 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/>.
\*---------------------------------------------------------------------------*/
#include "makeCombustionTypes.H"
#include "psiCombustionModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace combustionModels
{
makeCombustionTypes
(
infinitelyFastChemistry,
psiCombustionModel
);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

211
src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C
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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -58,10 +58,8 @@ void Foam::layerAdditionRemoval::checkDefinition()
{
if (!faceZoneID_.active())
{
FatalErrorIn
(
"void Foam::layerAdditionRemoval::checkDefinition()"
) << "Master face zone named " << faceZoneID_.name()
FatalErrorIn("void Foam::layerAdditionRemoval::checkDefinition()")
<< "Master face zone named " << faceZoneID_.name()
<< " cannot be found."
<< abort(FatalError);
}
@ -79,13 +77,15 @@ void Foam::layerAdditionRemoval::checkDefinition()
<< abort(FatalError);
}
if (topoChanger().mesh().faceZones()[faceZoneID_.index()].empty())
label nFaces = topoChanger().mesh().faceZones()[faceZoneID_.index()].size();
reduce(nFaces, sumOp<label>());
if (nFaces == 0)
{
FatalErrorIn
(
"void Foam::layerAdditionRemoval::checkDefinition()"
) << "Face extrusion zone contains no faces. "
<< " Please check your mesh definition."
FatalErrorIn("void Foam::layerAdditionRemoval::checkDefinition()")
<< "Face extrusion zone contains no faces. "
<< "Please check your mesh definition."
<< abort(FatalError);
}
@ -96,6 +96,7 @@ void Foam::layerAdditionRemoval::checkDefinition()
}
}
Foam::scalar Foam::layerAdditionRemoval::readOldThickness
(
const dictionary& dict
@ -115,21 +116,22 @@ void Foam::layerAdditionRemoval::clearAddressing() const
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::layerAdditionRemoval::layerAdditionRemoval
(
const word& name,
const label index,
const polyTopoChanger& mme,
const polyTopoChanger& ptc,
const word& zoneName,
const scalar minThickness,
const scalar maxThickness
const scalar maxThickness,
const Switch thicknessFromVolume
)
:
polyMeshModifier(name, index, mme, true),
faceZoneID_(zoneName, mme.mesh().faceZones()),
polyMeshModifier(name, index, ptc, true),
faceZoneID_(zoneName, ptc.mesh().faceZones()),
minLayerThickness_(minThickness),
maxLayerThickness_(maxThickness),
thicknessFromVolume_(thicknessFromVolume),
oldLayerThickness_(-1.0),
pointsPairingPtr_(NULL),
facesPairingPtr_(NULL),
@ -140,19 +142,22 @@ Foam::layerAdditionRemoval::layerAdditionRemoval
}
// Construct from dictionary
Foam::layerAdditionRemoval::layerAdditionRemoval
(
const word& name,
const dictionary& dict,
const label index,
const polyTopoChanger& mme
const polyTopoChanger& ptc
)
:
polyMeshModifier(name, index, mme, Switch(dict.lookup("active"))),
faceZoneID_(dict.lookup("faceZoneName"), mme.mesh().faceZones()),
polyMeshModifier(name, index, ptc, Switch(dict.lookup("active"))),
faceZoneID_(dict.lookup("faceZoneName"), ptc.mesh().faceZones()),
minLayerThickness_(readScalar(dict.lookup("minLayerThickness"))),
maxLayerThickness_(readScalar(dict.lookup("maxLayerThickness"))),
thicknessFromVolume_
(
dict.lookupOrDefault<Switch>("thicknessFromVolume", true)
),
oldLayerThickness_(readOldThickness(dict)),
pointsPairingPtr_(NULL),
facesPairingPtr_(NULL),
@ -189,11 +194,13 @@ bool Foam::layerAdditionRemoval::changeTopology() const
// Layer removal:
// When the min thickness falls below the threshold, trigger removal.
const faceZone& fz = topoChanger().mesh().faceZones()[faceZoneID_.index()];
const polyMesh& mesh = topoChanger().mesh();
const faceZone& fz = mesh.faceZones()[faceZoneID_.index()];
const labelList& mc = fz.masterCells();
const scalarField& V = topoChanger().mesh().cellVolumes();
const vectorField& S = topoChanger().mesh().faceAreas();
const scalarField& V = mesh.cellVolumes();
const vectorField& S = mesh.faceAreas();
if (min(V) < -VSMALL)
{
@ -206,63 +213,68 @@ bool Foam::layerAdditionRemoval::changeTopology() const
scalar avgDelta = 0;
scalar minDelta = GREAT;
scalar maxDelta = 0;
label nDelta = 0;
forAll(fz, faceI)
if (thicknessFromVolume_)
{
scalar curDelta = V[mc[faceI]]/mag(S[fz[faceI]]);
avgDelta += curDelta;
minDelta = min(minDelta, curDelta);
maxDelta = max(maxDelta, curDelta);
// Thickness calculated from cell volume/face area
forAll(fz, faceI)
{
scalar curDelta = V[mc[faceI]]/mag(S[fz[faceI]]);
avgDelta += curDelta;
minDelta = min(minDelta, curDelta);
maxDelta = max(maxDelta, curDelta);
}
nDelta = fz.size();
}
else
{
// Thickness calculated from edges on layer
const Map<label>& zoneMeshPointMap = fz().meshPointMap();
// Edges with only one point on zone
forAll(mc, faceI)
{
const cell& cFaces = mesh.cells()[mc[faceI]];
const edgeList cellEdges(cFaces.edges(mesh.faces()));
forAll(cellEdges, i)
{
const edge& e = cellEdges[i];
if (zoneMeshPointMap.found(e[0]))
{
if (!zoneMeshPointMap.found(e[1]))
{
scalar curDelta = e.mag(mesh.points());
avgDelta += curDelta;
nDelta++;
minDelta = min(minDelta, curDelta);
maxDelta = max(maxDelta, curDelta);
}
}
else
{
if (zoneMeshPointMap.found(e[1]))
{
scalar curDelta = e.mag(mesh.points());
avgDelta += curDelta;
nDelta++;
minDelta = min(minDelta, curDelta);
maxDelta = max(maxDelta, curDelta);
}
}
}
}
}
avgDelta /= fz.size();
reduce(minDelta, minOp<scalar>());
reduce(maxDelta, maxOp<scalar>());
reduce(avgDelta, sumOp<scalar>());
reduce(nDelta, sumOp<label>());
////MJ Alternative thickness determination
//{
// // Edges on layer.
// const Map<label>& zoneMeshPointMap = fz().meshPointMap();
//
// label nDelta = 0;
//
// // Edges with only one point on zone
// const polyMesh& mesh = topoChanger().mesh();
//
// forAll(mc, faceI)
// {
// const cell& cFaces = mesh.cells()[mc[faceI]];
// const edgeList cellEdges(cFaces.edges(mesh.faces()));
//
// forAll(cellEdges, i)
// {
// const edge& e = cellEdges[i];
//
// if (zoneMeshPointMap.found(e[0]))
// {
// if (!zoneMeshPointMap.found(e[1]))
// {
// scalar curDelta = e.mag(mesh.points());
// avgDelta += curDelta;
// nDelta++;
// minDelta = min(minDelta, curDelta);
// maxDelta = max(maxDelta, curDelta);
// }
// }
// else
// {
// if (zoneMeshPointMap.found(e[1]))
// {
// scalar curDelta = e.mag(mesh.points());
// avgDelta += curDelta;
// nDelta++;
// minDelta = min(minDelta, curDelta);
// maxDelta = max(maxDelta, curDelta);
// }
// }
// }
// }
//
// avgDelta /= nDelta;
//}
avgDelta /= nDelta;
if (debug)
{
@ -287,7 +299,6 @@ bool Foam::layerAdditionRemoval::changeTopology() const
}
// No topological changes allowed before first mesh motion
//
oldLayerThickness_ = avgDelta;
topologicalChange = false;
@ -315,7 +326,7 @@ bool Foam::layerAdditionRemoval::changeTopology() const
<< "Triggering layer removal" << endl;
}
triggerRemoval_ = topoChanger().mesh().time().timeIndex();
triggerRemoval_ = mesh.time().timeIndex();
// Old thickness looses meaning.
// Set it up to indicate layer removal
@ -347,7 +358,7 @@ bool Foam::layerAdditionRemoval::changeTopology() const
<< "Triggering layer addition" << endl;
}
triggerAddition_ = topoChanger().mesh().time().timeIndex();
triggerAddition_ = mesh.time().timeIndex();
// Old thickness looses meaning.
// Set it up to indicate layer removal
@ -377,9 +388,9 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
// Clear addressing. This also resets the addition/removal data
if (debug)
{
Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
<< " for object " << name() << " : "
<< "Clearing addressing after layer removal. " << endl;
Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
<< "for object " << name() << " : "
<< "Clearing addressing after layer removal" << endl;
}
triggerRemoval_ = -1;
@ -393,9 +404,9 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
// Clear addressing. This also resets the addition/removal data
if (debug)
{
Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
<< " for object " << name() << " : "
<< "Clearing addressing after layer addition. " << endl;
Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
<< "for object " << name() << " : "
<< "Clearing addressing after layer addition" << endl;
}
triggerAddition_ = -1;
@ -409,8 +420,8 @@ void Foam::layerAdditionRemoval::updateMesh(const mapPolyMesh&)
if (debug)
{
Pout<< "layerAdditionRemoval::updateMesh(const mapPolyMesh&) "
<< " for object " << name() << " : "
<< "Clearing addressing on external request. ";
<< "for object " << name() << " : "
<< "Clearing addressing on external request";
if (pointsPairingPtr_ || facesPairingPtr_)
{
@ -431,16 +442,14 @@ void Foam::layerAdditionRemoval::updateMesh(const mapPolyMesh&)
void Foam::layerAdditionRemoval::setMinLayerThickness(const scalar t) const
{
if
(
t < VSMALL
|| maxLayerThickness_ < t
)
if (t < VSMALL || maxLayerThickness_ < t)
{
FatalErrorIn
(
"void layerAdditionRemoval::setMinLayerThickness("
"const scalar t) const"
"void layerAdditionRemoval::setMinLayerThickness"
"("
"const scalar"
") const"
) << "Incorrect layer thickness definition."
<< abort(FatalError);
}
@ -451,15 +460,14 @@ void Foam::layerAdditionRemoval::setMinLayerThickness(const scalar t) const
void Foam::layerAdditionRemoval::setMaxLayerThickness(const scalar t) const
{
if
(
t < minLayerThickness_
)
if (t < minLayerThickness_)
{
FatalErrorIn
(
"void layerAdditionRemoval::setMaxLayerThickness("
"const scalar t) const"
"void layerAdditionRemoval::setMaxLayerThickness"
"("
"const scalar"
") const"
) << "Incorrect layer thickness definition."
<< abort(FatalError);
}
@ -475,7 +483,8 @@ void Foam::layerAdditionRemoval::write(Ostream& os) const
<< faceZoneID_ << nl
<< minLayerThickness_ << nl
<< oldLayerThickness_ << nl
<< maxLayerThickness_ << endl;
<< maxLayerThickness_ << nl
<< thicknessFromVolume_ << endl;
}
@ -490,6 +499,8 @@ void Foam::layerAdditionRemoval::writeDict(Ostream& os) const
<< token::END_STATEMENT << nl
<< " maxLayerThickness " << maxLayerThickness_
<< token::END_STATEMENT << nl
<< " thicknessFromVolume " << thicknessFromVolume_
<< token::END_STATEMENT << nl
<< " oldLayerThickness " << oldLayerThickness_
<< token::END_STATEMENT << nl
<< " active " << active()

16
src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.H
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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -65,6 +65,10 @@ class layerAdditionRemoval
//- Max thickness of extrusion layer. Triggers layer addition
mutable scalar maxLayerThickness_;
//- Switch to calculate thickness as volume/area
// If false, thickness calculated from edges
const bool thicknessFromVolume_;
//- Layer thickness from previous step
// Used to decide whether to add or remove layers
mutable scalar oldLayerThickness_;
@ -78,7 +82,7 @@ class layerAdditionRemoval
//- Layer removal trigger time index
mutable label triggerRemoval_;
//- Layer addition trigger time index
//- Layer addition trigger time index
mutable label triggerAddition_;
@ -120,6 +124,7 @@ class layerAdditionRemoval
//- Clear addressing
void clearAddressing() const;
// Helpers
//- Optionally read old thickness
@ -149,10 +154,11 @@ public:
(
const word& name,
const label index,
const polyTopoChanger& mme,
const polyTopoChanger& ptc,
const word& zoneName,
const scalar minThickness,
const scalar maxThickness
const scalar maxThickness,
const Switch thicknessFromVolume = true
);
//- Construct from dictionary
@ -161,7 +167,7 @@ public:
const word& name,
const dictionary& dict,
const label index,
const polyTopoChanger& mme
const polyTopoChanger& ptc
);

4
src/lagrangian/basic/particle/particle.C
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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -34,6 +34,8 @@ const Foam::scalar Foam::particle::trackingCorrectionTol = 1e-5;
const Foam::scalar Foam::particle::lambdaDistanceToleranceCoeff = 1e3*SMALL;
const Foam::scalar Foam::particle::minStepFractionTol = 1e5*SMALL;
namespace Foam
{
defineTypeNameAndDebug(particle, 0);

3
src/lagrangian/basic/particle/particle.H
View File

@ -307,6 +307,9 @@ public:
// for the denominator and numerator of lambda
static const scalar lambdaDistanceToleranceCoeff;
//- Minimum stepFraction tolerance
static const scalar minStepFractionTol;
// Constructors

6
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
View File

@ -652,7 +652,11 @@ void Foam::KinematicCloud<CloudType>::scaleSources()
template<class CloudType>
void Foam::KinematicCloud<CloudType>::preEvolve()
{
Info<< "\nSolving cloud " << this->name() << endl;
// force calculaion of mesh dimensions - needed for parallel runs
// with topology change due to lazy evaluation of valid mesh dimensions
label nGeometricD = mesh_.nGeometricD();
Info<< "\nSolving " << nGeometricD << "-D cloud " << this->name() << endl;
this->dispersion().cacheFields(true);
forces_.cacheFields(true);

4
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -110,11 +110,11 @@ void Foam::cloudSolution::read()
dict_.lookup("transient") >> transient_;
dict_.lookup("coupled") >> coupled_;
dict_.lookup("cellValueSourceCorrection") >> cellValueSourceCorrection_;
dict_.lookup("maxCo") >> maxCo_;
if (steadyState())
{
dict_.lookup("calcFrequency") >> calcFrequency_;
dict_.lookup("maxCo") >> maxCo_;
dict_.lookup("maxTrackTime") >> maxTrackTime_;
if (coupled_)

72
src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.C
View File

@ -69,13 +69,6 @@ bool Foam::CollidingParcel<ParcelType>::move
typename TrackData::cloudType::parcelType& p =
static_cast<typename TrackData::cloudType::parcelType&>(*this);
td.switchProcessor = false;
td.keepParticle = true;
const polyMesh& mesh = td.cloud().pMesh();
const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
const scalarField& V = mesh.cellVolumes();
switch (td.part())
{
case TrackData::tpVelocityHalfStep:
@ -92,70 +85,7 @@ bool Foam::CollidingParcel<ParcelType>::move
case TrackData::tpLinearTrack:
{
scalar tEnd = (1.0 - p.stepFraction())*trackTime;
const scalar dtMax = tEnd;
while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
{
// Apply correction to position for reduced-D cases
meshTools::constrainToMeshCentre(mesh, p.position());
const point start(p.position());
// Set the Lagrangian time-step
scalar dt = min(dtMax, tEnd);
// Remember which cell the parcel is in since this
// will change if a face is hit
const label cellI = p.cell();
const scalar magU = mag(p.U());
if (p.active() && magU > ROOTVSMALL)
{
const scalar d = dt*magU;
const scalar maxCo = td.cloud().solution().maxCo();
const scalar dCorr = min(d, maxCo*cbrt(V[cellI]));
dt *=
dCorr/d
*p.trackToFace(p.position() + dCorr*p.U()/magU, td);
}
tEnd -= dt;
p.stepFraction() = 1.0 - tEnd/trackTime;
// Avoid problems with extremely small timesteps
if (dt > ROOTVSMALL)
{
// Update cell based properties
p.setCellValues(td, dt, cellI);
if (td.cloud().solution().cellValueSourceCorrection())
{
p.cellValueSourceCorrection(td, dt, cellI);
}
p.calc(td, dt, cellI);
}
if (p.onBoundary() && td.keepParticle)
{
if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
{
td.switchProcessor = true;
}
}
p.age() += dt;
td.cloud().functions().postMove
(
p,
cellI,
dt,
start,
td.keepParticle
);
}
ParcelType::move(td, trackTime);
break;
}

21
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
View File

@ -271,6 +271,8 @@ bool Foam::KinematicParcel<ParcelType>::move
scalar tEnd = (1.0 - p.stepFraction())*trackTime;
const scalar dtMax = tEnd;
bool moving = true;
while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
{
// Apply correction to position for reduced-D cases
@ -281,12 +283,11 @@ bool Foam::KinematicParcel<ParcelType>::move
// Set the Lagrangian time-step
scalar dt = min(dtMax, tEnd);
// Remember which cell the parcel is in since this will change if
// a face is hit
// Cache the parcel current cell as this will change if a face is hit
const label cellI = p.cell();
const scalar magU = mag(U_);
if (p.active() && magU > ROOTVSMALL)
if (p.active() && moving && (magU > ROOTVSMALL))
{
const scalar d = dt*magU;
const scalar dCorr = min(d, maxCo*cbrt(V[cellI]));
@ -296,7 +297,19 @@ bool Foam::KinematicParcel<ParcelType>::move
}
tEnd -= dt;
p.stepFraction() = 1.0 - tEnd/trackTime;
scalar newStepFraction = 1.0 - tEnd/trackTime;
if
(
mag(p.stepFraction() - newStepFraction)
< particle::minStepFractionTol
)
{
moving = false;
}
p.stepFraction() = newStepFraction;
// Avoid problems with extremely small timesteps
if (dt > ROOTVSMALL)

11
src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -57,7 +57,7 @@ Foam::StandardWallInteraction<CloudType>::StandardWallInteraction
"StandardWallInteraction<CloudType>::StandardWallInteraction"
"("
"const dictionary&, "
"CloudType& cloud"
"CloudType&"
")"
) << "Unknown interaction result type "
<< interactionTypeName
@ -174,10 +174,11 @@ bool Foam::StandardWallInteraction<CloudType>::correct
(
"bool StandardWallInteraction<CloudType>::correct"
"("
"typename CloudType::parcelType&, "
"const polyPatch&, "
"const label, "
"bool&, "
"vector&"
"bool& keepParticle, "
"const scalar, "
"const tetIndices&"
") const"
) << "Unknown interaction type "
<< this->interactionTypeToWord(interactionType_)

75
src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C
View File

@ -36,7 +36,7 @@ Foam::PilchErdman<CloudType>::PilchErdman
:
BreakupModel<CloudType>(dict, owner, typeName),
B1_(0.375),
B2_(0.236)
B2_(0.2274)
{
if (!this->defaultCoeffs(true))
{
@ -90,58 +90,79 @@ bool Foam::PilchErdman<CloudType>::update
scalar& massChild
)
{
scalar semiMass = nParticle*pow3(d);
scalar We = 0.5*rhoc*sqr(Urmag)*d/sigma;
// Weber number - eq (1)
scalar We = rhoc*sqr(Urmag)*d/sigma;
// Ohnesorge number - eq (2)
scalar Oh = mu/sqrt(rho*d*sigma);
scalar Wec = 6.0*(1.0 + 1.077*pow(Oh, 1.6));
// Critical Weber number - eq (5)
scalar Wec = 12.0*(1.0 + 1.077*pow(Oh, 1.6));
if (We > Wec)
{
// We > 1335, wave crest stripping
// We > 2670, wave crest stripping - eq (12)
scalar taubBar = 5.5;
if (We < 1335)
if (We < 2670)
{
if (We > 175.0)
if (We > 351)
{
// sheet stripping
taubBar = 0.766*pow(2.0*We - 12.0, 0.25);
// sheet stripping - eq (11)
taubBar = 0.766*pow(We - 12.0, 0.25);
}
else if (We > 22.0)
else if (We > 45)
{
// Bag-and-stamen breakup
taubBar = 14.1*pow(2.0*We - 12.0, -0.25);
// bag-and-stamen breakup - eq (10)
taubBar = 14.1*pow(We - 12.0, 0.25);
}
else if (We > 9.0)
else if (We > 18)
{
// Bag breakup
taubBar = 2.45*pow(2.0*We - 12.0, 0.25);
// bag breakup - eq (9)
taubBar = 2.45*pow(We - 12.0, 0.25);
}
else if (We > 6.0)
else if (We > 12)
{
// Vibrational breakup
taubBar = 6.0*pow(2.0*We - 12.0, -0.25);
// vibrational breakup - eq (8)
taubBar = 6.0*pow(We - 12.0, -0.25);
}
else
{
// no break-up
taubBar = GREAT;
}
}
scalar rho12 = sqrt(rhoc/rho);
scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*taubBar*taubBar);
// velocity of fragmenting drop - eq (20)
scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*sqr(taubBar));
// maximum stable diameter - eq (33)
scalar Vd1 = sqr(1.0 - Vd/Urmag);
Vd1 = max(Vd1, SMALL);
scalar Ds = 2.0*Wec*sigma/(Vd1*rhoc*sqr(Urmag));
scalar A = Urmag*rho12/d;
scalar dStable = Wec*sigma/(Vd1*rhoc*sqr(Urmag));
scalar taub = taubBar/A;
if (d < dStable)
{
// droplet diameter already stable = no break-up
// - do not update d and nParticle
return false;
}
else
{
scalar semiMass = nParticle*pow3(d);
scalar frac = dt/taub;
// invert eq (3) to create a dimensional break-up time
scalar taub = taubBar*d/(Urmag*rho12);
// update the droplet diameter according to the rate eq. (implicitly)
d = (d + frac*Ds)/(1.0 + frac);
// update droplet diameter according to the rate eq (implicitly)
scalar frac = dt/taub;
d = (d + frac*dStable)/(1.0 + frac);
// correct the number of particles to conserve mass
nParticle = semiMass/pow3(d);
// correct the number of particles to conserve mass
nParticle = semiMass/pow3(d);
}
}
return false;

22
src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.H
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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ Class
Foam::PilchErdman
Description
secondary breakup model
Particle secondary breakup model, based on the reference:
@verbatim
Pilch, M. and Erdman, C.A.
@ -35,6 +35,24 @@ Description
Int. J. Multiphase Flows 13 (1987), 741-757
@endverbatim
The droplet fragment velocity is described by the equation:
\f[
V_d = V sqrt(epsilon)(B1 T + B2 T^2)
\f]
Where:
V_d : fragment velocity
V : magnitude of the relative velocity
epsilon : density ratio (rho_carrier/rho_droplet)
T : characteristic break-up time
B1, B2 : model input coefficients
The authors suggest that:
compressible flow : B1 = 0.75*1.0; B2 = 3*0.116
incompressible flow : B1 = 0.75*0.5; B2 = 3*0.0758
\*---------------------------------------------------------------------------*/
#ifndef PilchErdman_H

15
src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C
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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -49,7 +49,7 @@ const ThermoType& Foam::multiComponentMixture<ThermoType>::constructSpeciesData
template<class ThermoType>
void Foam::multiComponentMixture<ThermoType>::correctMassFractions()
{
// It changes Yt patches to "calculated"
// Multiplication by 1.0 changes Yt patches to "calculated"
volScalarField Yt("Yt", 1.0*Y_[0]);
for (label n=1; n<Y_.size(); n++)
@ -57,6 +57,17 @@ void Foam::multiComponentMixture<ThermoType>::correctMassFractions()
Yt += Y_[n];
}
if (mag(max(Yt).value()) < ROOTVSMALL)
{
FatalErrorIn
(
"void Foam::multiComponentMixture<ThermoType>::"
"correctMassFractions()"
)
<< "Sum of mass fractions is zero for species " << this->species()
<< exit(FatalError);
}
forAll(Y_, n)
{
Y_[n] /= Yt;