Files
openfoam/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
2012-03-30 08:53:47 +01:00

214 lines
5.4 KiB
C

/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "ParticleErosion.H"
// * * * * * * * * * * * * * Protectd Member Functions * * * * * * * * * * * //
template<class CloudType>
Foam::label Foam::ParticleErosion<CloudType>::applyToPatch
(
const label globalPatchI
) const
{
forAll(patchIDs_, i)
{
if (patchIDs_[i] == globalPatchI)
{
return i;
}
}
return -1;
}
template<class CloudType>
void Foam::ParticleErosion<CloudType>::write()
{
if (QPtr_.valid())
{
QPtr_->write();
}
else
{
FatalErrorIn("void Foam::ParticleErosion<CloudType>::write()")
<< "QPtr not valid" << abort(FatalError);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class CloudType>
Foam::ParticleErosion<CloudType>::ParticleErosion
(
const dictionary& dict,
CloudType& owner
)
:
CloudFunctionObject<CloudType>(dict, owner, typeName),
QPtr_(NULL),
patchIDs_(),
p_(readScalar(this->coeffDict().lookup("p"))),
psi_(this->coeffDict().template lookupOrDefault<scalar>("psi", 2.0)),
K_(this->coeffDict().template lookupOrDefault<scalar>("K", 2.0))
{
const wordList allPatchNames = owner.mesh().boundaryMesh().names();
wordList patchName(this->coeffDict().lookup("patches"));
labelHashSet uniquePatchIDs;
forAllReverse(patchName, i)
{
labelList patchIDs = findStrings(patchName[i], allPatchNames);
if (patchIDs.empty())
{
WarningIn
(
"Foam::ParticleErosion<CloudType>::ParticleErosion"
"("
"const dictionary&, "
"CloudType& "
")"
) << "Cannot find any patch names matching " << patchName[i]
<< endl;
}
uniquePatchIDs.insert(patchIDs);
}
patchIDs_ = uniquePatchIDs.toc();
// trigger ther creation of the Q field
preEvolve();
}
template<class CloudType>
Foam::ParticleErosion<CloudType>::ParticleErosion
(
const ParticleErosion<CloudType>& pe
)
:
CloudFunctionObject<CloudType>(pe),
QPtr_(NULL),
patchIDs_(pe.patchIDs_),
p_(pe.p_),
psi_(pe.psi_),
K_(pe.K_)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class CloudType>
Foam::ParticleErosion<CloudType>::~ParticleErosion()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CloudType>
void Foam::ParticleErosion<CloudType>::preEvolve()
{
if (QPtr_.valid())
{
QPtr_->internalField() = 0.0;
}
else
{
const fvMesh& mesh = this->owner().mesh();
QPtr_.reset
(
new volScalarField
(
IOobject
(
this->owner().name() + "Q",
mesh.time().timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("zero", dimVolume, 0.0)
)
);
}
}
template<class CloudType>
void Foam::ParticleErosion<CloudType>::postPatch
(
const parcelType& p,
const label patchI,
const label patchFaceI
)
{
const label localPatchI = applyToPatch(patchI);
if (localPatchI != -1)
{
const fvMesh& mesh = this->owner().mesh();
// patch-normal direction
vector nw = p.currentTetIndices().faceTri(mesh).normal();
// particle direction of travel
const vector& U = p.U();
// quick reject if particle travelling away from the patch
if ((-nw & U) < 0)
{
return;
}
nw /= mag(nw);
const scalar magU = mag(U);
const vector Udir = U/magU;
// determine impact angle, alpha
const scalar alpha = mathematical::pi/2.0 - acos(nw & Udir);
const scalar coeff = p.nParticle()*p.mass()*sqr(magU)/(p_*psi_*K_);
scalar& Q = QPtr_->boundaryField()[patchI][patchFaceI];
if (tan(alpha) < K_/6.0)
{
Q += coeff*(sin(2.0*alpha) - 6.0/K_*sqr(sin(alpha)));
}
else
{
Q += coeff*(K_*sqr(cos(alpha))/6.0);
}
}
}
// ************************************************************************* //