ENH: Adding per-patch wall interaction.

Rearranged data storage in existing wall model.  Improved stability.
This commit is contained in:
graham
2010-05-25 19:46:58 +01:00
parent bcdc6291ca
commit f4fd4c00dc
11 changed files with 604 additions and 54 deletions

View File

@ -39,7 +39,7 @@ inline bool Foam::WallCollisionRecord<Type>::match
// Using the new data as the acceptance criterion
scalar cosAcceptanceAngle = magpRel/radius;
if (cosAcceptanceAngle > 1.0)
if (cosAcceptanceAngle > 1.0 + SMALL)
{
Info<< "pRel_ " << pRel_ << " " << magpRel_ << nl
<< "pRel " << pRel << " " << magpRel << nl

View File

@ -36,6 +36,7 @@ License
#include "PairSpringSliderDashpot.H"
#include "WallSpringSliderDashpot.H"
#include "WallLocalSpringSliderDashpot.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -73,6 +74,13 @@ License
WallSpringSliderDashpot, \
KinematicCloud, \
ParcelType \
); \
\
makeWallModelType \
( \
WallLocalSpringSliderDashpot, \
KinematicCloud, \
ParcelType \
);

View File

@ -37,6 +37,7 @@ License
#include "PairSpringSliderDashpot.H"
#include "WallSpringSliderDashpot.H"
#include "WallLocalSpringSliderDashpot.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -99,9 +100,16 @@ License
KinematicCloud, \
ParcelType, \
ThermoType \
); \
\
makeWallModelThermoType \
( \
WallLocalSpringSliderDashpot, \
KinematicCloud, \
ParcelType, \
ThermoType \
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

View File

@ -100,11 +100,11 @@ Foam::PairSpringSliderDashpot<CloudType>::PairSpringSliderDashpot
)
),
volumeFactor_(1.0),
useEquivalentSize_(Switch(this->dict().lookup("useEquivalentSize")))
useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
{
if (useEquivalentSize_)
{
volumeFactor_ = readScalar(this->dict().lookup("volumeFactor"));
volumeFactor_ = readScalar(this->coeffDict().lookup("volumeFactor"));
}
scalar nu = this->owner().constProps().poissonsRatio();

View File

@ -58,9 +58,6 @@ class PairSpringSliderDashpot
// the same Poisson's ratio and Young's modulus
scalar Gstar_;
//- Poisson's ratio of both particles
scalar sigma_;
//- alpha-coefficient, related to coefficient of restitution
scalar alpha_;
@ -93,7 +90,7 @@ class PairSpringSliderDashpot
scalar volumeFactor_;
//- Switch to control use of equivalent size particles. Used
// because cbrt calculation is very expensive.
// because the calculation can be very expensive.
bool useEquivalentSize_;

View File

@ -0,0 +1,354 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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/>.
\*---------------------------------------------------------------------------*/
#include "WallLocalSpringSliderDashpot.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template <class CloudType>
void Foam::WallLocalSpringSliderDashpot<CloudType>::findMinMaxProperties
(
scalar& rMin,
scalar& rhoMax,
scalar& UMagMax
) const
{
rMin = VGREAT;
rhoMax = -VGREAT;
UMagMax = -VGREAT;
forAllConstIter(typename CloudType, this->owner(), iter)
{
const typename CloudType::parcelType& p = iter();
// Finding minimum diameter to avoid excessive arithmetic
scalar dEff = p.d();
if (useEquivalentSize_)
{
dEff *= cbrt(p.nParticle()*volumeFactor_);
}
rMin = min(dEff, rMin);
rhoMax = max(p.rho(), rhoMax);
UMagMax = max
(
mag(p.U()) + mag(p.omega())*dEff/2,
UMagMax
);
}
// Transform the minimum diameter into minimum radius
// rMin = dMin/2
rMin /= 2.0;
}
template <class CloudType>
void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
(
typename CloudType::parcelType& p,
const point& site,
const WallSiteData<vector>& data,
scalar pREff
) const
{
// wall patch index
label wPI = patchMap_[data.patchIndex()];
// data for this patch
scalar Estar = Estar_[wPI];
scalar Gstar = Gstar_[wPI];
scalar alpha = alpha_[wPI];
scalar b = b_[wPI];
scalar mu = mu_[wPI];
vector r_PW = p.position() - site;
vector U_PW = p.U() - data.wallData();
scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);
vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;
scalar etaN = alpha*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
vector fN_PW =
rHat_PW
*(kN*pow(normalOverlapMag, b) - etaN*(U_PW & rHat_PW));
p.f() += fN_PW;
vector USlip_PW =
U_PW - (U_PW & rHat_PW)*rHat_PW
+ (p.omega() ^ (pREff*-rHat_PW));
scalar deltaT = this->owner().mesh().time().deltaTValue();
vector& tangentialOverlap_PW =
p.collisionRecords().matchWallRecord(-r_PW, pREff).collisionData();
tangentialOverlap_PW += USlip_PW*deltaT;
scalar tangentialOverlapMag = mag(tangentialOverlap_PW);
if (tangentialOverlapMag > VSMALL)
{
scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar;
scalar etaT = etaN;
// Tangential force
vector fT_PW;
if (kT*tangentialOverlapMag > mu*mag(fN_PW))
{
// Tangential force greater than sliding friction,
// particle slips
fT_PW = -mu*mag(fN_PW)*USlip_PW/mag(USlip_PW);
tangentialOverlap_PW = vector::zero;
}
else
{
fT_PW =
-kT*tangentialOverlapMag
*tangentialOverlap_PW/tangentialOverlapMag
- etaT*USlip_PW;
}
p.f() += fT_PW;
p.torque() += (pREff*-rHat_PW) ^ fT_PW;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template <class CloudType>
Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
(
const dictionary& dict,
CloudType& cloud
)
:
WallModel<CloudType>(dict, cloud, typeName),
Estar_(),
Gstar_(),
alpha_(),
b_(),
mu_(),
patchMap_(),
maxEstarIndex_(-1),
collisionResolutionSteps_
(
readScalar
(
this->coeffDict().lookup("collisionResolutionSteps")
)
),
volumeFactor_(1.0),
useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
{
if (useEquivalentSize_)
{
volumeFactor_ = readScalar(this->coeffDict().lookup("volumeFactor"));
}
scalar pNu = this->owner().constProps().poissonsRatio();
scalar pE = this->owner().constProps().youngsModulus();
const polyMesh& mesh = cloud.mesh();
const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
patchMap_.setSize(bMesh.size(), -1);
DynamicList<label> wallPatchIndices;
forAll(bMesh, patchI)
{
if (isA<wallPolyPatch>(bMesh[patchI]))
{
wallPatchIndices.append(bMesh[patchI].index());
}
}
label nWallPatches = wallPatchIndices.size();
Estar_.setSize(nWallPatches);
Gstar_.setSize(nWallPatches);
alpha_.setSize(nWallPatches);
b_.setSize(nWallPatches);
mu_.setSize(nWallPatches);
scalar maxEstar = -GREAT;
forAll(wallPatchIndices, wPI)
{
const dictionary& patchCoeffDict
(
this->coeffDict().subDict(bMesh[wallPatchIndices[wPI]].name())
);
patchMap_[wallPatchIndices[wPI]] = wPI;
scalar nu = dimensionedScalar
(
patchCoeffDict.lookup("poissonsRatio")
).value();
scalar E = dimensionedScalar
(
patchCoeffDict.lookup("youngsModulus")
).value();
Estar_[wPI] = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);
Gstar_[wPI] = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));
alpha_[wPI] =
dimensionedScalar(patchCoeffDict.lookup("alpha")).value();
b_[wPI] = dimensionedScalar(patchCoeffDict.lookup("b")).value();
mu_[wPI] = dimensionedScalar(patchCoeffDict.lookup("mu")).value();
if (Estar_[wPI] > maxEstar)
{
maxEstarIndex_ = wPI;
maxEstar = Estar_[wPI];
}
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template <class CloudType>
Foam::WallLocalSpringSliderDashpot<CloudType>::~WallLocalSpringSliderDashpot()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CloudType>
Foam::scalar Foam::WallLocalSpringSliderDashpot<CloudType>::pREff
(
const typename CloudType::parcelType& p
) const
{
if (useEquivalentSize_)
{
return p.d()/2*cbrt(p.nParticle()*volumeFactor_);
}
else
{
return p.d()/2;
}
}
template<class CloudType>
bool Foam::WallLocalSpringSliderDashpot<CloudType>::controlsTimestep() const
{
return true;
}
template<class CloudType>
Foam::label Foam::WallLocalSpringSliderDashpot<CloudType>::nSubCycles() const
{
if (!(this->owner().size()))
{
return 1;
}
scalar rMin;
scalar rhoMax;
scalar UMagMax;
findMinMaxProperties(rMin, rhoMax, UMagMax);
// Note: pi^(7/5)*(5/4)^(2/5) = 5.429675
scalar minCollisionDeltaT =
5.429675
*rMin
*pow(rhoMax/(Estar_[maxEstarIndex_]*sqrt(UMagMax) + VSMALL), 0.4)
/collisionResolutionSteps_;
return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
}
template<class CloudType>
void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
(
typename CloudType::parcelType& p,
const List<point>& flatSitePoints,
const List<WallSiteData<vector> >& flatSiteData,
const List<point>& sharpSitePoints,
const List<WallSiteData<vector> >& sharpSiteData
) const
{
scalar pREff = this->pREff(p);
forAll(flatSitePoints, siteI)
{
evaluateWall
(
p,
flatSitePoints[siteI],
flatSiteData[siteI],
pREff
);
}
forAll(sharpSitePoints, siteI)
{
// Treating sharp sites like flat sites
evaluateWall
(
p,
sharpSitePoints[siteI],
sharpSiteData[siteI],
pREff
);
}
}
// ************************************************************************* //

View File

@ -0,0 +1,184 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 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::WallLocalSpringSliderDashpot
Description
Forces between particles and walls, interacting with a spring,
slider, damper model
\*---------------------------------------------------------------------------*/
#ifndef WallLocalSpringSliderDashpot_H
#define WallLocalSpringSliderDashpot_H
#include "WallModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class WallLocalSpringSliderDashpot Declaration
\*---------------------------------------------------------------------------*/
template<class CloudType>
class WallLocalSpringSliderDashpot
:
public WallModel<CloudType>
{
// Private data
//- Effective Young's modulus value
scalarList Estar_;
//- Effective shear modulus value
scalarList Gstar_;
//- alpha-coefficient, related to coefficient of restitution
scalarList alpha_;
//- Spring power (b = 1 for linear, b = 3/2 for Hertzian)
scalarList b_;
//- Coefficient of friction in for tangential sliding
scalarList mu_;
//- Mapping the patch index to the model data
labelList patchMap_;
//- Index of the maximum value of Estar_
label maxEstarIndex_;
//- The number of steps over which to resolve the minimum
// harmonic approximation of the collision period
scalar collisionResolutionSteps_;
//- Volume factor for determining the equivalent size of a
// parcel where nParticles is not 1. The equivalent size of
// the parcel is
// parcelEquivVolume = volumeFactor*nParticles*p.volume()
// so
// parcelEquivD = cbrt(volumeFactor*nParticles)*p.d()
// + When volumeFactor = 1, the particles are compressed
// together so that the equivalent volume of the parcel is
// the sum of the constituent particles
// + When volumeFactor = 3*sqrt(2)/pi, the particles are
// close packed, but uncompressed.
// + When volumeFactor > 3*sqrt(2)/pi, the particles loosely
// grouped.
// 3*sqrt(2)/pi = 1.350474 is the volume factor for close
// packing, i.e pi/(3*sqrt(2)) is the maximum close packing
// factor
scalar volumeFactor_;
//- Switch to control use of equivalent size particles. Used
// because the calculation can be very expensive.
bool useEquivalentSize_;
// Private Member Functions
//- Find the appropriate properties for determining the minimum
//- allowable timestep
void findMinMaxProperties
(
scalar& rMin,
scalar& rhoMax,
scalar& vMagMax
) const;
//- Calculate the wall interaction for a parcel at a given site
void evaluateWall
(
typename CloudType::parcelType& p,
const point& site,
const WallSiteData<vector>& data,
scalar pREff
) const;
public:
//- Runtime type information
TypeName("WallLocalSpringSliderDashpot");
// Constructors
//- Construct from dictionary
WallLocalSpringSliderDashpot(const dictionary& dict, CloudType& cloud);
//- Destructor
virtual ~WallLocalSpringSliderDashpot();
// Member Functions
//- Return the volumeFactor
inline scalar volumeFactor() const
{
return volumeFactor_;
}
//- Return the effective radius for a particle for the model
virtual scalar pREff(const typename CloudType::parcelType& p) const;
//- Whether the WallModel has a timestep limit that will
// require subCycling
virtual bool controlsTimestep() const;
//- For WallModels that control the timestep, calculate the
// number of subCycles needed to satisfy the minimum
// allowable timestep
virtual label nSubCycles() const;
//- Calculate the wall interaction for a parcel
virtual void evaluateWall
(
typename CloudType::parcelType& p,
const List<point>& flatSitePoints,
const List<WallSiteData<vector> >& flatSiteData,
const List<point>& sharpSitePoints,
const List<WallSiteData<vector> >& sharpSiteData
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "WallLocalSpringSliderDashpot.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -58,6 +58,14 @@ Foam::WallModel<CloudType>::owner() const
}
template<class CloudType>
CloudType&
Foam::WallModel<CloudType>::owner()
{
return owner_;
}
template<class CloudType>
const Foam::dictionary& Foam::WallModel<CloudType>::dict() const
{

View File

@ -111,6 +111,9 @@ public:
//- Return the owner cloud object
const CloudType& owner() const;
//- Return non-const access to the owner cloud object
CloudType& owner();
//- Return the dictionary
const dictionary& dict() const;

View File

@ -76,19 +76,15 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
typename CloudType::parcelType& p,
const point& site,
const WallSiteData<vector>& data,
scalar pNu,
scalar pE,
scalar pREff,
scalar Estar,
scalar kN,
scalar Gstar
scalar kN
) const
{
vector r_PW = p.position() - site;
vector U_PW = p.U() - data.wallData();
scalar normalOverlapMag = pREff - mag(r_PW);
scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);
vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
@ -115,7 +111,7 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
if (tangentialOverlapMag > VSMALL)
{
scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar;
scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar_;
scalar etaT = etaN;
@ -156,8 +152,8 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
)
:
WallModel<CloudType>(dict, cloud, typeName),
E_(dimensionedScalar(this->coeffDict().lookup("youngsModulus")).value()),
nu_(dimensionedScalar(this->coeffDict().lookup("poissonsRatio")).value()),
Estar_(),
Gstar_(),
alpha_(dimensionedScalar(this->coeffDict().lookup("alpha")).value()),
b_(dimensionedScalar(this->coeffDict().lookup("b")).value()),
mu_(dimensionedScalar(this->coeffDict().lookup("mu")).value()),
@ -169,12 +165,30 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
)
),
volumeFactor_(1.0),
useEquivalentSize_(Switch(this->dict().lookup("useEquivalentSize")))
useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
{
if (useEquivalentSize_)
{
volumeFactor_ = readScalar(this->dict().lookup("volumeFactor"));
volumeFactor_ = readScalar(this->coeffDict().lookup("volumeFactor"));
}
scalar nu = dimensionedScalar
(
this->coeffDict().lookup("poissonsRatio")
).value();
scalar E = dimensionedScalar
(
this->coeffDict().lookup("youngsModulus")
).value();
scalar pNu = this->owner().constProps().poissonsRatio();
scalar pE = this->owner().constProps().youngsModulus();
Estar_ = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);
Gstar_ = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));
}
@ -225,17 +239,11 @@ Foam::label Foam::WallSpringSliderDashpot<CloudType>::nSubCycles() const
findMinMaxProperties(rMin, rhoMax, UMagMax);
scalar pNu = this->owner().constProps().poissonsRatio();
scalar pE = this->owner().constProps().youngsModulus();
scalar Estar = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu_))/E_);
// Note: pi^(7/5)*(5/4)^(2/5) = 5.429675
scalar minCollisionDeltaT =
5.429675
*rMin
*pow(rhoMax/(Estar*sqrt(UMagMax) + VSMALL), 0.4)
*pow(rhoMax/(Estar_*sqrt(UMagMax) + VSMALL), 0.4)
/collisionResolutionSteps_;
return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
@ -252,17 +260,9 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
const List<WallSiteData<vector> >& sharpSiteData
) const
{
scalar pNu = this->owner().constProps().poissonsRatio();
scalar pE = this->owner().constProps().youngsModulus();
scalar pREff = this->pREff(p);
scalar Estar = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu_))/E_);
scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;
scalar GStar = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu_ - sqr(nu_))/E_));
scalar kN = (4.0/3.0)*sqrt(pREff)*Estar_;
forAll(flatSitePoints, siteI)
{
@ -271,12 +271,8 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
p,
flatSitePoints[siteI],
flatSiteData[siteI],
pNu,
pE,
pREff,
Estar,
kN,
GStar
kN
);
}
@ -289,12 +285,8 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
p,
sharpSitePoints[siteI],
sharpSiteData[siteI],
pNu,
pE,
pREff,
Estar,
kN,
GStar
kN
);
}
}

View File

@ -50,11 +50,11 @@ class WallSpringSliderDashpot
{
// Private data
//- Young's modulus of the wall
scalar E_;
//- Effective Young's modulus value
scalar Estar_;
//- Poisson's ratio of the wall
scalar nu_;
//- Effective shear modulus value
scalar Gstar_;
//- alpha-coefficient, related to coefficient of restitution
scalar alpha_;
@ -109,12 +109,8 @@ class WallSpringSliderDashpot
typename CloudType::parcelType& p,
const point& site,
const WallSiteData<vector>& data,
scalar pNu,
scalar pE,
scalar pREff,
scalar Estar,
scalar kN,
scalar Gstar
scalar kN
) const;